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Abstract 


In this paper ; we propose a new technique for the numerical treat- 
ment of external flow problems with oscillatory behavior of the solu- 
tion in time. Specifically, we consider the case of unbounded 
compressible viscous plane flow past a finite body (airfoil). Oscilla- 
tions of the flow in time may be caused by the time -periodic injection of 
fluid into the boundary layer y which in accordance with experimental 
data , may essentially increase the performance of the airfoil To con- 
duct the actual computations , we have to somehow restrict the original 
unbounded domain, that is, to introduce an artificial ( external ) bound- 
ary and to further consider only a finite computational domain. Conse- 
quently, we will need to formulate some artificial boundary conditions 
(ABC's) at the introduced external boundary . The ABC's we are aim- 
ing to obtain must meet a fundamental requirement. One should be 
able to uniquely complement the solution calculated inside the finite 
computational domain to its infinite exterior so that the original prob- 
lem is solved within the desired accuracy. Our construction of such 
ABC's for oscillating flows is based on an essential assumption: the 
Navier-Stokes equations can be linearized in the far field against the 
free-stream background. To actually compute the ABC's, we represent 
the far-field solution as a Fourier series in time and then apply the Dif- 
ference Potentials Method (DPM) ofV. S. Ryaben'kii. This paper con- 
tains a general theoretical description of the algorithm for setting the 
DPM-based ABC's for time-periodic external flows. Based on our 
experience in implementing analogous ABC's for steady-state prob- 
lems (a simpler case), we expect that these boundary conditions will 
become an effective tool for constructing robust numerical methods to 
calculate oscillatory flows. 

1. Introduction 

The numerical study of problems originally formulated on unbounded domains requires the imple- 
mentation of special techniques for the “treatment of infinity” (which is necessitated by the restricted 
facilities of modem computers). One of the corresponding techniques is based on an artificial truncation 
of the original infinite domain, which implies that one must set special boundary conditions at the exter- 
nal (artificial) boundary of the newly formed finite computational domain. The aim of this paper is to 
describe the theoretical foundations for constructing such artificial boundary conditions (ABC’s) for the 
computation of certain unsteady external flows. 

Before proceeding to the actual description of the problem, let us first define the concept of exact 
ABC's. Namely, exact ABC’s are the boundary conditions that enable one to uniquely complement the 
solution of the “truncated problem” to the unbounded exterior of the computational domain so that the 
original problem is solved. The exact ABC’s usually appear to be nonlocal for steady-state problems in 
space and for time-dependent problems in both space and time. 

Let us emphasize that our main objective in this paper is to construct special boundary conditions 
that would model (and in the ideal case equivalently replace) the exterior part of the problem, i.e., the 
part we eliminate by truncation. Many examples of such boundary conditions can be found in compre- 
hensive reviews by Givoli. (See refs. 1 and 2.) This formulation differs from another well-known prob- 
lem related to setting the boundary conditions for numerical algorithms, namely, to construct such 
boundary conditions that would ensure well-posedness of the truncated problem and stability of the 



integration process in time. In fact, these two formulations are not completely independent. For exam- 
ple, the issue of well-posedness for certain classes of (local) ABC’s was thoroughly investigated by 
Gustafsson. (See refs. 3-5.) On the other hand, a group of very delicate questions related to the issue of 
long-time stability is studied by Carpenter, Gottlieb, and Abarbanel in reference 6 (for some specific 
boundary-value problems). The issue of connections between the (highly accurate nonlocal) boundary 
conditions that “model the infinity” and the boundary conditions that ensure the long-time stability will 
be an interesting subject for a future investigation. 

In this paper, we consider an unbounded compressible viscous flow past a finite body or configura- 
tion of bodies (e.g., single-element or multi-element airfoil). The behavior of the flow in time is 
assumed to be oscillatory. We must emphasize that while talking about the oscillatory time behavior we 
mean that some alternating (time-periodic) influence is exerted on the flow (e.g., see experimental work 
by Seifert, et al. in ref. 7) and expect that those frequencies that are connected to this influence will 
dominate in the solution. We expect that this circumstance will enable us to construct the ABC’s with- 
out taking into account any other time-dependent effects. From a mathematical standpoint, this case fills 
an intermediate position between the steady-state and true unsteady flows. 

The steady-state case is relatively simple compared with time-dependent flows. In reference 8, we 
have constructed the ABC’s for calculating external viscous compressible steady-state flows. These 
boundary conditions were based on the concept of far-field linearization and on the application of the 
Difference Potentials Method (DPM) of Ryaben’kii. (See refs. 9 and 10.) The ABC’s (ref. 8) differ only 
slightly from the exact ABC’s (a more rigorous formulation of the latter statement may be found in 
ref. 8); therefore, the ABC’s (ref. 8) turn out to be global in space. However, practical implementation 
of these boundary conditions is fairly easy. (See refs. 1 1 and 12.) They were used along with the Navier- 
Stokes code by Jameson, Schmidt, Turkel, and Swanson (refs. 13-15) for computing different external 
flows. Numerical experiments show that the global DPM-based ABC’s (ref. 8) provide high accuracy of 
computations, as well as fast convergence of the multigrid iteration procedure to a steady state. (See 
refs. 1 1 and 12.) The computational cost of boundary conditions (refs. 8, 1 1, and 12) is not high in com- 
parison with the total cost of the original procedure. (See refs. 13-15.) Generally, the numerical algo- 
rithm we used for integrating the Navier-Stokes equations became more robust (in comparison with the 
standard procedure (refs. 13-15)) if supplemented by the DPM-based ABC’s. (See ref. 8.) 

Additionally, we would like to emphasize that the ABC’s (ref. 8) were constructed specially for the 
steady-state problem and on the basis of stationary governing equations, independent of any specific 
technique for solving the stationary equations inside the computational domain. In practical computa- 
tions, we use multigrid iterations (refs. 13-15) for calculating the steady-state solutions in references 1 1 
and 12. In doing so, we set the ABC’s (ref. 8) on each iteration on the upper time level. Of course, the 
boundary data on the intermediate stage of the iteration procedure (i.e., until we achieve the steady 
state) are not necessarily consistent with the formal “stationary” treatment of the far field. However, 
treating the “time-intermediate” boundary data as if it were already steady has been found effective in 
computational practice. (See refs. 1 1 and 12.) We are going to use a similar idea for the time-periodic 
case studied in this paper. 

True unsteady flows are much more complicated in terms of both theoretical analysis and practical 
calculations. In general, the exact ABC’s for unsteady problems will be nonlocal in both space and time. 
Therefore, the corresponding computational cost may appear to be rather high. This is also true for the 
global DPM-based boundary conditions which can be constructed as close to the exact ones as desired. 
The corresponding general theory for unsteady problems is contained in work by Ryaben’kii. (See 
ref. 16.) 

However, an intermediate case of oscillatory time behavior must be less expensive in terms of 
required computer resources since the global character of the ABC’s in time will obviously be restricted 
by the value of one period. Moreover, the theoretical analysis of this case based on the usage of 
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the Fourier representation in time also appears to be less complicated than the general one from 
reference 16 since in our analysis we actually reduce the time-dependent problem to a family of steady- 
state problems. 

On the other hand, don’t assume that the oscillating flow is a particular and, therefore, an unimpor- 
tant case. For example, experiments (ref. 7) show that the time-periodic injection of fluid into the turbu- 
lent boundary layer may increase its resistance to adverse pressure gradients without separation. This 
implies an essential improvement of airfoil performance, up to 60 percent for high (post stall) angles of 
attack, according to reference 7. The phenomenon was observed on different geometries (original 
NACA0015 airfoil, the same airfoil with the deflected flap, and some others), which leads us to believe 
that it may be effectively used in aircraft design. Therefore, an accurate numerical investigation of the 
phenomenon becomes an important issue, and an accurate procedure for setting the ABC’s must be one 
of the principle elements of any computational algorithm used for such an investigation. 

The previous example is probably not a unique one where the time-periodic treatment of flow in the 
far field might be relevant. In general, for the oscillatory case we propose the following construction of 
ABC’s. First, linearize the governing equations in the far field. Then, assuming that the time period is 
initially prescribed, apply the Fourier transform in time and obtain a family of steady-state problems 
(where the unknowns are amplitudes). The latter problems are then treated by means of the DPM. (See 
refs. 9 and 10.) The central idea of the DPM-based approach is to equivalently replace the problem for- 
mulated on a domain by a certain operator equation formulated on its boundary. For each one of the 
foregoing steady-state problems (note, the family of these problems is parameterized by the frequency, 
i.e., by the dual Fourier variable), this replacement results in an operator equation formulated at the 
artificial boundary of the computational domain (i.e., connecting the boundary values of the solution). 
The operator involved (a projection) is somewhat analogous to the boundary pseudodifferential opera- 
tors introduced by Calderon. (See ref. 1 7.) Because of the equivalence to the exterior linear problem, the 
previously-mentioned operator equation (more precisely, the entire family of these equations) can be 
considered a desirable exact ABC (limited only by the accuracy of far-field linearization) for the prob- 
lem solved inside the computational domain. In other words, this operator equation adequately takes 
into account the structure of the solution from outside the computational domain, which might also be 
called the exact transfer of boundary conditions from infinity. (See ref. 16.) 

We actually develop the DPM-based ABC’s for the already discrete formulation of the problem. In 
doing so, the set of the frequencies involved is obviously finite. Therefore, we can actually compute the 
corresponding boundary operator for each one of the steady-state problems arising after the Fourier 
transform in time. Then, for reasons of numerical convenience, we represent the solution to the linear- 
ized exterior problem in the form of generalized potential. (See refs. 9 and 10.) The density of general- 
ized potential serves as an unknown function in the previously-mentioned operator equation. By using 
the generalized potential to set the ABC’s we gain more generality from a geometric standpoint. More- 
over, we can easily match the solutions of the interior nonlinear problem and the exterior linear problem 
when conducting practical computations. (We need to actually calculate the generalized potential only 
in some neighborhood of the computational domain, to be discussed later.) Finally, applying the inverse 
Fourier transform, we obtain ABC’s in a matrix form, which enables easy practical implementation. In 
fact, the entire procedure may be thought of as solving the linearized problem outside the computational 
domain and then using the obtained solution to close the “truncated system” that is solved inside the 
computational domain. 

To conclude this introduction, let us point out an analogy to the previously investigated steady-state 
case. (See ref. 8.) Namely, we are looking here for a solution to the unsteady problem that is defined on 
an initially prescribed time interval (one period) and that meets the periodicity condition in time (at least 
in the far field). To develop the ABC’s for this case, we solve a certain linear problem in the far field (by 
means of the DPM). The latter problem is also formulated for the time interval of one period. The 
ABC’s for the time-periodic case are basically constructed independent of any specific technique for 
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integrating the Navier-Stokes equations inside the computational domain (as the ABC’s (ref. 8) were 
constructed irrespective of any specific way for actual computation of the steady state). Based on the 
assumption of periodicity in time, these boundary conditions simply close the system that is solved 
inside the computational domain; the closure is constructed for the time interval of one period. In prac- 
tice, however, achieving the true oscillatory regime may require long-time computational runs that 
cover many periods. During this long-time integration, each moment we need to update the external 
boundary data using the ABC’s (i.e., each time step) we treat the flow as it were already time-periodic 
(in some generalized sense, see section 3). In so doing, the boundary conditions should guarantee only 
the desirable far-field behavior of the solution. This behavior is actually determined by the condition 
that all perturbations vanish at infinity (as in refs. 8, 11, and 12 when we were treating the external 
boundary data on each iteration as already steady and demanding that the ABC’s ensure the decrease of 
the solution to the linearized problem at infinity). 

This paper is organized as follows. In section 2, we describe the basic formulations of the problems. 
Specifically, in subsection 2.1 we describe a geometric setup typical for the numerical solution of exter- 
nal flow problems, i.e., configurations of the finite computational domain and its infinite exterior. In 
this subsection, we also introduce the flow equations (parabolized Navier-Stokes) and linearize them in 
the far field against the constant free-stream background. In so doing, we obtain a coupled problem, 
which is nonlinear inside the finite computational domain and linear outside it. Then, assuming that the 
period of oscillating motion is known, we Fourier transform the exterior linear system with respect to 
time and obtain an equivalent family of steady-state systems. These steady-state systems must be solved 
as a part of the solution to the aforementioned coupled problem. However, we do not solve them 
directly since the corresponding domain is still infinite. Instead, we equivalently replace each of the 
exterior linear steady-state systems by the generalized Calderon pseudodifferential equation formulated 
at the external boundary of the computational domain. To calculate the pseudodifferential operation 
(projection) we need a special auxiliary problem that is first formulated on the entire plane for the lin- 
earized thin-layer equations (after the Fourier transform in time) with a certain compactly supported 
right-hand side. Solvability of this auxiliary problem (in the sense of tempered distributions) is studied 
in subsection 2.2. Then, in subsection 2.3, we show how one can replace the original auxiliary problem 
formulated on an unbounded domain (entire plane) by a new problem formulated on some rectangle so 
that the solutions of the two problems are in a certain sense close to each other. 

Section 3 of this paper is devoted to numerics. In subsection 3.1, we introduce a finite-difference 
scheme that approximates the linearized thin-layer equations. Since we discretize the equations not only 
in space but also in time, we now get a finite (discrete) series instead of the original infinite Fourier 
series which implies that the family of steady-state systems to be solved outside the computational 
domain becomes finite as well. In subsection 3.2, we construct a difference analogue to the auxiliary 
problem on the rectangle, describe the numerical algorithm for its solution (referring to our previous 
work for some details) and briefly address our somewhat non-standard concept of convergence for the 
solutions of the difference auxiliary problem. Finally, in subsection 3.3 we show how one uses the 
recently formulated difference auxiliary problem and obtains difference analogue to the Calderon 
boundary pseudodifferential projection. Using this difference boundary projection and also calculating 
the generalized difference potential, we actually compute the nonlocal DPM-based ABC’s. The ABC’s 
are first obtained in the Fourier variables and then, after implementing the inverse transform, in the 
physical variables as well. Finally, section 4 contains some conclusions and possible generalizations. 

2. Basic Formulations 

2.1. Governing Equations and Geometric Setup 

Let us start with the parabolized Navier-Stokes equations, which are the same as the thin-layer 
equations for two dimensions (see ref. 18 by Anderson, Tannehill, and Pletcher): 
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3p 3p^ 3pv 
3/ 3;t 3y 


= 0 


du du du dp 1 3 3 u 
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P di + pU di + P % + Ty ~ Re 3 


'a/ K “a* 


fdu + _ 1 

fdu'f 4 fdu) 2 y d del 

V3x 3 y) Re 

+ 3H^J + Prd-y%\ 


( 1 ) 


Here, x and v denote the Cartesian coordinates, u and v denote the Cartesian velocity projections, p 
denotes the density, p denotes the pressure, £ denotes the internal energy, p denotes the viscosity, and y 
denotes the ratio of specific heats. To derive the last of equations (1), we assume that the gas is perfect 
and that the Prandtl number Pr = pc^/K is constant (k is the heat conduction coefficient). We denote the 
free-stream parameters, specifically, u 0 , v 0 , p 0 , p 0 , £ 0 , and jOq, by the subscript “0.” We additionally 
assume that Vq = 0 and w 0 > 0, which does not imply any loss of generality. The system (1) is written in 
dimensionless form. The following scales were used to obtain dimensionless quantities: w 0 was used for 

velocity; p 0 for density, Pou 0 2 for pressure, u 0 2 for internal energy, Pq for viscosity, characteristic size L 
(typically, airfoil chord) for all distances, and L/uq for time. The factor MRe that multiplies the viscous 


terms in equations (1) arises from the nondimensionalization. Here, Re - 


Po“o L 

Po 


number. 


is the Reynolds 


Note that in our previous work (refs. 8, 1 1, and 12) we used the full Navier-Stokes equations to con- 
struct the ABC’s for steady-state problems. In this paper, we are going to use the thin-layer system 
(eqs. (1)). This system appears to apply quite well to the description of certain viscous flows (ref. 18), in 
particular, the far-field flows that we are studying hereafter. Moreover, for the thin-layer system 
(eqs. (1)) we can justify some results on the solvability of its linearized counterpart on /? 2 , which is 
important for the general justification of our construction of ABC’s. Finally, the usage of equations (1) 
instead of the full Navier-Stokes equations may save an appreciable amount of computer resources, as 
will be seen from further consideration. 


Let us now assume that the actual values of w, v, /?, p, e, and (X in the far field only slightly deviate 
from the corresponding free-stream parameters. For dimensionless quantities, this means 


P = 1 +P 


U = 1 + U V = V |i - 1 + ji 

e = [(y- l + e 


p = (yM$)~' + p 


( 2 ) 


where 


p « 1 u « 1 v « 1 p « 1 p « (yM^)~ l 
e«[(y-l)yM2]-> 

Here, M 0 = u 0 (yp 0 /p 0 )~ l/2 is the Mach number at infinity, which is always assumed to be less than 
unity. By substituting expressions (2) into equations (1) and retaining only the first-order terms with 
respect to small perturbations u , v , p , p , e , and ji , we obtain the following system of linear partial 
differential equations with constant coefficients: 
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0 
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0 

0 
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Re 

0 4/3 

0 

0 
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0 0 

yPr ~ 1 

-Pt~ 1 Mq 2 


(3b) 


The system (3) is the linearization of equations (1) against the free-stream background. We omit the 
tilde in equations (3) since we are going to deal only with linear equations in perturbations henceforth. 

Additionally, we used the equation of state e = — - (more precisely, its linearization 


e = 


Y-l 


7^o 


p ) to eliminate internal energy from equations (3). 


We have mentioned that equations (3) will be used for the description of fluid motion in the far 
field. Let us now define a general geometric setup for the problem under consideration. The original 
Navier-Stokes equations are integrated on a grid (e.g., C-type) generated around the airfoil; this grid 
covers the finite computational domain which is denoted D in hereafter. (See fig. 1.) We henceforth 
assume that the linearization (eqs. (3)) is valid outside the computational domain D in , i.e., on its com- 
plement D ex . (See fig. 1 .) This assumption is true for large computational domains, i.e., far enough from 
the immersed body. As we approach the airfoil, the possibility of linearization in D ex can always be ver- 
ified a posteriori by analyzing the corresponding computational results (as was done in refs. 1 1 and 12 
for the steady-state case). 

To integrate the Navier-Stokes equations on the grid inside D in , we use some finite-dimensional 
approximation of these equations. The actual type of the resulting discrete operator (i.e., finite- 
difference, finite-element, etc.) is not that important from the standpoint of constructing the ABC’s; for 
definiteness we assume that the Navier-Stokes equations are integrated by means of a finite-difference 
scheme. To begin with, we also suppose that this scheme is fully explicit in time. We may think that we 
already know the solution for the time level i on the entire grid, in particular, / = 0 implies the initial 
data. When we advance one time step, i.e., calculate the solution for the level t M by means of the 
scheme, we cannot obtain this solution for the whole grid since some nodes located near the external 
boundary of D in will be missing. The actual location of missing nodes depends on the specific structure 
of the scheme stencil. For example, a typical central-difference second-order approximation to the spa- 
tial part of the Navier-Stokes operator on a structured grid requires a 3 X 3 stencil. Using such a spatial 
approximation combined with an explicit integration procedure in time, we can obtain the solution on 
the level at all nodes, except for those that belong to the outermost coordinate row of the grid (des- 
ignated Tj in fig. 1). To advance the next time step (t l+2 ) we will have to somehow determine these 
missing values of the solution on the level This will be done by means of solving the linearized sys- 
tem in D ex (i.e., by representing its solution in the form of the generalized potential for each Fourier 
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Figure 1. Configuration of domains. 

mode). In other words, using the solution to equations (3) in D ex , we close the system of difference 
equations inside the computational domain D in . The closure we obtain is actually the desirable ABC’s. 

Note that in the case of steady-state problems the ABC’s (ref. 8) were also used to close the subdef- 
inite system of difference equations inside D in . As previously mentioned, boundary conditions (ref. 8) 
were implemented in references 1 1 and 12 together with the pseudo-time iteration procedure for achiev- 
ing the steady state. (See refs. 13-15.) From an algorithmic standpoint, this approach is almost the same 
as the true integration in time, so the ABC’s (ref. 8) were applied on the upper time level for each itera- 
tion. However, there is an essential difference between the approach in reference 8 and the technique to 
be described in this paper. Namely, the former is intended only for the treatment of steady-state prob- 
lems and is based on the linearized stationary equations, and the latter will take into account the previ- 
ous evolution of the solution in time. 

Additionally, let us note that in the case of implicit schemes we also need ABC’s that will complete 
the system of difference equations inside D in . Indeed, while integrating the Navier-Stokes equations by 
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means of an implicit scheme one has to solve a certain discrete system on the upper time level (/ +1 ), 
whereas the data from the lower time level(s) play the role of forcing terms. This system will obviously 
be subdefinite unless we specify additional relations that connect the values of unknowns in the grid 
nodes located near the external boundary. In particular, for the previously-mentioned example of a 
structured grid and central differences on the 3 x 3 spatial stencil, these additional relations (i.e., the 
ABC’s) should connect the values of the solution at the penultimate (the curve T in fig. 1 ) and outermost 
rows of grid nodes. (See also refs. 8, 1 1, and 12.) Including the missing relations provided by the ABC’s 
into the system solved on the upper time level, we close this system and then advance the next time step. 

Let us now provide an exact formulation of the problem. First, we select those nodes of the grid 
where the solution can no longer be determined by the scheme but must be obtained by means of special 
additional relations (i.e., by means of the ABC’s). We designate this set of nodes v } . Second, we select 
those nodes of the grid where we need to know the solution in order to obtain it on Vj with the help of 
the ABC’s. The latter set is designated v. Both v and Vj will depend on the structure of the specific sten- 
cil. In particular, for the 3x3 stencil on a structured grid, v and V| correspond to the penultimate and 
outermost rows of grid nodes, respectively. (Also see refs. 8, 1 1 , and 12.) Without loss of generality, we 
assume that the artificial boundary T (see fig. 1) is formed by the penultimate row of nodes v, so that all 
nodes Vj that form the curve Y j (see fig. 1) already belong to D ex (i.e., to the “linear zone”). 

Then, we designate the time period by T '. Clearly, we can further consider our problem for the time 
interval [0,7] without loss of generality. We will also need the following brief notations: 
DJ x = n ex X [0, T] , Dj n = D in x [0, T) , r T = r x [0, T] , and T[ = Tj x [0, T] . The closure of the 
finite-difference system in Df which we are looking for and which should be provided by the ABC’s, 
is actually a set of relations expressing u| r r in terms of some data specified on Y T . As previously men- 
tioned, these relations will be based on the solution to the linearized system (3) in Dj x . The latter sys- 
tem is supplemented (on Dj x ) by the periodicity condition in time, 

u Uo = u l, = 7- i(x,y)zD ex ) (4) 

and the free-stream condition at infinity, 

u — > 0 as x 2 + y 2 — > oo (t e [0,T]) (5) 

The choice of the data on T^that “drive” the ABC’s is closely connected to the concept of clear trace , 
delineated in references 9 and 10. The question of the possible proper constructions of clear traces for 
equations (3) may require a special thorough investigation in addition to the general analysis from refer- 
ences 9 and 10; such an investigation is not a direct subject of this paper. Therefore, we will not com- 
ment on this question in our further discussion, we only point out the actual construction we use. 
Namely, let us first represent the vector function u(x,y,r) in the form of a Fourier series in time for any 
space point ft,y), 


n = 2n 

m ^ A - l fit „ 

U(*, y, t) = 2L u (x,y)e 

n = — 


( 6 ) 


where 


«"(*, y) = t)e 


. 2n 
-ini— 

T dt 


(n = 0, ±1,±2, ...) 


(7) 
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Instead of considering equations (3a), (4), and (5) Dj x , we henceforth consider D ex the family of “sta- 
tionary” systems. 


/(O n Ca n + D^ + F^ + H^ = 0 (n = 0, ±1, ±2, ...) (8) 

parameterized by the frequency cd„ = 2nn/T, « = 0, ±1,±2, and supplemented at infinity by the 
boundary conditions 


u"Uy) ->0 as x 2 + y 2 — > oo ( n = 0,±1,±2, ...) (9) 

which directly follow from formula (5). The matrices C, D, F, and H in the system (8) are the same as in 
formula (3b). 


For each frequency co n we consider the pair of functions 




Up, 


(specified on T) as the data that 


“drive” the ABC’s; here, ^ is the normal to T. (Note that if the interior solution is already computed by 

du n 

means of the scheme inside Dj , then up and can be easily calculated.) 

in i ell, 


Our ultimate goal will be to provide a full classification of those and only those functions 


*n 0 \ p 

v r> 


V 




(defined on V) that generate a solution u"(jt, y) to system (8) (with boundary conditions (9) defined on 
D ex and such that its trace on T coincides with the “source” function itself, i.e.. 


au n> 


( ^ 

a« dv r 

l K) 

r 

v ^ ) 


(10) 


-n 

Vr ’^ 


can be 


As will be seen from further consideration, the corresponding set of functions 
described as an image of a certain boundary projection operator. In other words, the functions 
will satisfy some boundary operator equation with projection. (The equation of this type was 


f ^ n\ 

~ n dvp 

Vr * ac 


mentioned in the introduction as the one equivalent to the linearized exterior problem.) Let us designate 
the corresponding projection operator by Pp (we actually construct this operator in section 3). Then, 


specifying any function 


ac 


J r> 


(from inside D in ), we apply Pp and consider the projection 


P« 


p ac , 


a n dvp 

Vp 


V 


’ as 


as the right-hand side in equality (10) for the problem (eqs. (8)-( 10)). 


After solving the problem (eqs. (8)— (10)) on D ex , we find the trace of its solution on r t (i.e., on Vj), 
which in turn enables us to obtain the missing boundary relations that close the system of difference 
equations inside Df n . These relations (i.e., the ABC’s) are derived using the inverse Fourier transform. 
They can be symbolically written as 
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(re [0, T]) 


( 11 ) 


u = 



p 

( 

V 





— 

P" oP«oR 

u", ^ 



ex r 

l v a? J 



where the operator R represents some (smooth) interpolation of the discrete functions along the curve T, 
and the operator P* x involves the calculation of the generalized potential to solve the problem 
(eqs. (8)-(10)). The specific structure of all operators from formula (11) will be delineated in section 3, 
where we actually construct their discrete counterparts. 


Let us make a few important remarks. First, to formally close the system solved in Dj n , we have to 

obtain additional relations between the values of the unknowns on P^and on T[. Such relations would 
provide ABC’s that are completely independent of any specific numerical procedure employed inside 
Dj . However, to simplify our task and at the same time only slightly compromise the previously- 

mentioned independence, we take into account that we almost always integrate the Navier-Stokes equa- 
tions step-by-step in time (explicitly or implicitly). Therefore, we do not have to construct such ABC’s 
that would connect the values of the solution at v and at Vj for all time moments t e [0, T ] . It suffices to 
determine w, v, /?, and p at Vj only for t-T (i.e., at the upper time level) since for all previous moments 
these values have been determined when calculating previous time steps. Moreover, the formulation of 
the problem (eqs. (8)-(10)), where the right-hand side from equality (10) belongs to the projection 


image. 


> \ 

A n dvp 
v r 


V 




e Im Pp , assumes that these data are a result of operating by on the Fourier trans- 


form 


u" 

p a; 


of some time-periodic function. However, in conducting the step-by-step integration in 


time, the actual data 




may not be periodic until we achieve a true oscillatory regime. There- 


fore, as mentioned in the introduction, any time we use the ABC’s we implement a certain generalized 
treatment of the external flow as being already time-periodic. Namely, instead of the true boundary data 


3u r 


u 


p 


at T 1 , we use the best approximation of this data by periodic functions in the sense of least 
squares. This approach will be delineated in section 3, which is devoted to numerics. 


Second, we are unable to directly solve the problem (eqs. (8)-(10)) on D ex since the domain is infi- 
nite. Handling of this problem will require the additional truncation. Recall that we have already trun- 
cated the original infinite domain and have obtained D in \ now we also truncate D ex in order to get a new 
linear problem formulated on a finite domain, and therefore, available for solution on the computer. 
This issue is addressed in subsection 2.3. 


Third, we certainly will not solve the problem (eqs. (8)— (10)) every time we need to obtain a closed 
system inside D in (i.e., each time step). Instead, using the linearity of the problem, we will specify some 
basis in the space of boundary data and solve the problem (eqs. (8)-(10)) one time for each basis func- 
tion. This approach will enable us to obtain the ABC’s in matrix form, which is very convenient for 
practical computing. (Also see refs. 8, 1 1, and 12.) 

Ultimately we will deal only with the finite-difference formulations and, consequently, with the 
finite Fourier series (instead of the infinite series (6), see section 3). In so doing, the discretization in 
time for the linearized exterior problem in Dj x should not necessarily coincide with the one used for the 
Navier-Stokes scheme inside Dj n . A more convenient method may be to use interpolation in time, 
which was previously proposed in reference 16. 


10 



Finally, let us mention that since we need to know the solution on T for the whole period T to 
restore the solution on Vj, the first few time steps (until the total time reaches T) will require some spe- 
cial treatment. It might be based on the usage of either a larger grid or some other external boundary 
conditions for the initial stage of integration in time. 

We now proceed to the actual construction of the operators involved in formula (1 1 ). This construc- 
tion will be essentially the same for all wavenumbers n (n is contained as a parameter in the correspond- 
ing expressions hereafter). 

As was mentioned before, the computation of the ABC’s (eq. (11)) consists of two stages. First we 
apply the projection Pjt to provide the proper boundary data (right-hand side of equality (10) for the 
problem (eqs. (8)— (1 0)). Then we find the solution to the problem (eqs. (8)— (10)) in the form of the gen- 
eralized potential (operator P * x ). Both of these stages will require the application of the DPM. (See 
refs. 9 and 10.) In particular, it appears that the computation of Pfl and P£ x requires the solution of the 
same auxiliary problem ( AP ) described in sections 2.2 and 2.3 for the continuous formulation and in 
section 3.2 for the difference formulation. This AP is actually the main element of the DPM-based 
approach. The Green operator of the AP plays in the theory of generalized potentials approximately the 
same role as the Green function (or the fundamental solution) plays in classical potential theory. (See 
refs. 9 and 10.) The AP is formulated on the entire plane (;t,y) for the inhomogeneous counterpart of 
system (8) with a certain compactly supported right-hand side f" = (/". /2> /"> /") (to be specified 
later on). Namely, we will need to solve the following system, 


" dx dy dy 2 


= r 


( 12 ) 


on R 2 , suppf"(x. y ) c D jn , and we will require that the solution be unique in the class of functions van- 
ishing at infinity. In other words, system (12) is supplemented by the following boundary condition. 


u"(*, y)^0 as x 2 + y 2 oo (13) 

which is the same as boundary conditions (9). 

Once we are able to solve the AP (eqs. (12) and (13)), we can construct the boundary operator P{1 , 
properly formulate the problem (eqs. (8)— (10)), and finally obtain its solution in the form of a general- 
ized potential. This is actually a very brief description of our DPM-based approach; it will be delineated 
in section 3 for the discrete formulation of the problem. Now we will investigate the solvability of the 
AP (eqs. (12) and (13)). 


2.2. Solvability of Linearized Problem on Entire Plane 

We will look for the solution to the AP (eqs. (12) and (13)) in the space of tempered distributions 
g' (see ref. 19 by Hormander or ref. 20 by Vladimirov), which is a conjugate space to the space g of all 
infinitely smooth functions defined on R 2 that decrease at infinity with all their derivatives faster than 
any power of ( x 2 + y 2 )“ 1/2 . Take the Fourier transform 


a | 

u"(£,ll) = —j j u n (x,y)e- i ^-^ydxdy 


(14a) 


^ ,T1) = 2^j jr(x,y)e- i ^- i ^dxdy 


(14b) 
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of both sides of system (12) and represent the result in the form of a matrix equation, 


Qu = f 


( 15 ) 


Note that in system (15) and henceforth in this subsection we drop the superscript n to simplify the nota- 
tions. Then, the symbol Q (eq. (15)) is given by 


Q = j'toC + iE,D + i'qF -ti 2 H 

it, m 

n 2 

0 


i«o + 5) + 2- 


0 

0 




0 

«•$ 

1 


1(0) + ^) 
0 

0 


RePrMl 


(16) 


We first show that system (15) is solvable in Q' . For the time being, we do not need any restrictive 
assumptions in regard to f ; as previously mentioned, f is compactly supported (suppf c: D in ), and with- 
out loss of generality we may think that f is absolutely integrable on R 2 (f e L l (R 2 )). Then, its Fourier 

transform f is bounded and continuous on R ; consequently, if we formally write down the solution to 
system (15) as 


u = Q _I f (17) 

then the properties of the right-hand side in equality (17) are fully determined by the inverse symbol 
Q _1 . Indeed, it is well known (ref. 20) that if the right-hand side of equality (17) is locally absolutely 
integrable on R then it defines the tempered distribution, i.e., the generalized function from g\ The lat- 
ter will coincide (in the sense of distributions) with the classical function Q _I (^> r|)f (^, r\) everywhere 
on /? 2 , except for the set of singularities of Q _1 (^, r|)f (^, r\) (if any). Since in our case the function 
f (^, T|) is continuous and bounded on /? 2 , then it suffices to determine whether the function Q -1 (£„ q) 
belongs to 4oc(* 2 )- 

To do this, we have to find all singularities of Q _1 (^, T|). Calculating the determinant of Q(£, r|), 
we obtain 


Q&r\) = 


_ ((0 + %) A + ± &(? + ^ 2 ) + (a± + l JL)_ 1 _^H 4 

1 M 2 ^ ^ Re 2 U 3 Pr J 3 M 2 Pn 


M$PrRe 2 


MgPrRe 2 


i ,{ (*> + 4) W f7 . Y 
Re 13 Pr 


(to + 4)^ 2 q 2 f4 1 

M^Re \3 Pr 


ft ,1 1 (CQ + ^nV . n 4 
b Pr) M 2 Re l Pr) 3 


4 (M + ^)T) 6 
1 PrRe 3 


(18) 


Here C, and T| are the variables and w, M 0 , Re, Pr, and y are parameters. We emphasize that both vari- 
ables \ and rj are supposed to be real (see formulas (14)); however, the coefficients of Q( rj) are, gen- 
erally speaking, complex. Thus, to find singular points of the symbol Q (eq. (16)), one has to find the 
real roots of (J(C, T|) (eq. (18)), which actually implies to find common real roots of two polynomials, 
Q(%, r|) and 3<2(J;, T|)- First, the point (£ = -co, r| = 0) is clearly one of such common roots. Then, we 
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note that 3<2(^, r|) turns into zero on the two entire straight lines, ^ = -(0 and r\ = 0. Moreover, 
"H) has no other roots on the line £, = -to, except for T) = 0. Further, substituting r\ = 0 into the 
equation T|) = 0 (see formula (18)) and assuming that £, * -0), we find the following two roots of 

^ , C0A/ n 

^l) = 0 that belong to the line r\ = 0: £,1 = t — — and = — . We also observe that if 

1 M Q 1 " 4 " M Q 


CO = 0 (which corresponds to the steady-state flows), then all three roots, (-CD, 0), 


f coM 0 

y**o 


,0 


\ 


j 


, and 


f 

V 


coM n ^ 

— 0 

1 + Af 0 ’ J 


, merge into one. 


In an attempt to find other real roots (if any) of Q(^r\) (eq. (18)), we divide the equation 
3(?(^, T|) = 0 by (cd + ^)r ] 2 /Re. (This is possible since we have already proven that no other zeros exist 
on the two lines £, = -to and r\ = 0, except for those already found.) The resulting equation, 


(cd + 



4 rr 
-y — l - — 
3 PrRe 2 


= 0 


(19) 


is of fourth order, and taking into account that the equation 9*g(^, r\) = 0 (see formula (18)) is of sixth 
order, we conclude that the polynomial Q ( £,, r\) may have not more than a finite number of isolated real 
roots in total (three of which have already been found). We emphasize here that this property (finite 
number of isolated real roots) presents an essential difference between the problem under investigation 
and classical acoustics problems in which the viscous terms in the governing equations are usually 
neglected. Namely, for the acoustics equations (i.e., linearized Euler equations) the singular points of 
the symbol are no longer isolated. They usually form a curve on the plane R 2 which may cause notice- 
able difficulties with justification of the uniqueness of solution. These difficulties are similar to those 
that arise in studying the Helmholtz equation, which may be referred to as describing acoustics in the 
stationary medium. We do not deal with Helmholtz-like equations in this paper; we only note that con- 
trary to the acoustics case, system (12) is presumably easier from this standpoint since the proof of 
uniqueness appears to be elementary. (See proposition 4.) 

Since equation (19) is of second order with respect to ^ we can resolve it for each T| and obtain 
explicit function(s) £,3 = ^(q). Because we are interested only in real solutions, we have to consider a 
few different cases. 


First, assume that co ^ 0. Then, rewrite equation (19) as 


s 2 


z + x. 

3 Pr 




■V 

PrRe 2 


= 0 


(20) 


and observe that, if = ^ + j Q + ~ ) , then equation (20) degenerates and therefore has a 

unique real solution ^ 0) = ^ 0) (q) for any T\. If M l > (| + ^-) Q + ■%- ) ' , then we can easily make 
sure that the discriminant 


D = 4® 2 


(MY- 

V3 Pr) 


Z + X__L(4 + ± 

l 3 Pr M 2 V3 Pr 


\3 Pr) mA Pr) 3' Pri 


PrRe 2 


( 21 ) 
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is always positive, which means that equation (19) has two different real solutions, £,3^ = ^^(rj) and 

= £®(T]) for any T|. If ^ + ~p~ r ) ’ ^ en conc *ition D>0 (see formula (21)) 

imposes certain restrictions on T|. Namely, we have 



where 

D = 1 fi , if 16 T” 2 (1 I 1 Y 4 I 1 f 7 I 1 1 f 4 I 1 V 

1 M^V Pr) 3 PrRe 2 M$\3 PrA3 PrAJ> Pr Pr) 

Therefore, in this case the real solutions to equation (19), £3 * = ^,3 \t) ) and ^ 2) = ^ 2) (q), exist only 
forr| within the above range. (See inequality (22)). 

Now consider the case co = 0 (which corresponds to the steady-state problem). From equation (19), 
we easily derive 


S 3 Pr 

__l (t+±: 

\ _ 

ni(i + r 

M^V Pr y 

I+ 4 y ^ 


(23) 

A/^\3 Pr, 

) ~ 

> VPrRe 2 


Equation (23) has real 

solutions, 


= tghn) 

and ^ 2) 

II 

kT ft 

N 

-3 

only for 


Ml > (j + /r)(j + p”) ' Otherwise, we conclude that the equation 3 {?(£>, T|) = 0 for co = 0 has no 

other real roots except for (J; = 0, r\ = 0) and therefore, the same is true for the equation £)(£;, r\) = 0. 


In practice, we have calculated explicit symbolic expressions for the functions 

^3 ^ = ^>3 ) ( r l)’ an d ^3^ = ^3 2 ) (ti) using Mathematica. (See ref. 21 by Wolfram.) (These expressions 
are not presented here because they are fairly cumbersome.) Then, substituting the functions 
^3°^ = ^3 0) (r|), £3? = £3 ^(T|), and = ^^(T]) into the second equation, T]) = 0, we obtain 

the algebraic equations with respect to only one variable t|. Clearly the above equations (which are dif- 
ferent for the different solutions, £,3°^ = ^3 °\t|), £ 3^ = ^3 \r\), and £3^ = ^^(r))) may have real 
root(s) if and only if the original equation Q(^, rj) = 0 has other real zero(s) besides those that have 

( a>M 0 ] ( (oM 0 ) 

already been found, (-co, 0), - — — , 0 , and -- , 0 . Therefore, we finally have reduced the 

1 - M n 1 + A/ n 

v u J V u / 

question about the real zeros of the equation Q(^ y r|) = 0 to the question about the real root(s) of certain 
algebraic equations of one variable. 

Regrettably, the resulting equations (after the substitution of £3^ = ^3 0) (r|), £3* = ^3 ^(rj ) , and 
= £^(r|) into T l) = 0) appear to be too complicated for obtaining general expressions for 

their real root(s). However, we may implement the following semi -numerical approach which provides 
fairly convincing results. 

First, note that the case co = 0 seems to be the simplest one. This case actually admits rigorous anal- 
ysis without doing any simplifying assumptions. As previously mentioned, equation (23) has no real 
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solutions for Afg < + p^j (which implies that the determinant (eq. (18)) has no real roots); 

for Af o ~ (j + 7>~rj<3> + Pr) ec l uat ‘ on ( 23 ) degenerates and any pair (£,, q) of the kind £, is arbitrary, 
q = 0 is its root. Substituting this root into q) = 0 (see eq. (18)), we obtain ^ 3 4 ( 1 /M\ - 1 ) = 0, 

which yields ^ = 0. Therefore, we did not find any new real zero. For Mq > , equa- 

tion (23) has two different real solutions for any q; moreover, ^ = ^ 0 (q) = -^ 2) = -^ 2) (q). 
Since all powers of c, in q) are even, we do not need to separately consider = £3 *(q) and 

^3* = Substituting ^g 1 ' 2 * = ^3’ 2) (q) into 9?0(^, q) = 0, we obtain the following eighth- 

order equation with respect to q: aq 8 + br] 6 + cq 4 = 0, where the coefficients a , b, and c are obviously 
real. The explicit expressions for a , b y and c were obtained by means of Mathematica (see ref. 21); we 
do not present them here because they are cumbersome. However, using these expressions we can prove 
that a > 0, b > 0, and c > 0. Then, it becomes clear that there are no other real roots except for the one 
we have already found, r\ = 0 (which also yields £, = 0). Indeed, the equation <rr | 4 + br \ 2 + c = 0 for 
a> 0, b > 0, and c > 0 may have only essentially complex roots r|. Therefore, we conclude that for 
0) = 0 the symbol Q (eq. (16)) has only one singular point (^ = 0, T| = 0). 

Recall that all equations under study generally depend on five real parameters, co, Mq, Re , Pr , and y. 
To simplify our task, we fix the values of some of these parameters. Let us set y = 1.4 (two-atom gas) 
and Pr - 0.72 (air). This choice of values for the ratio of specific heats and for the Prandtl number, 
respectively, is most frequently used since it is closely related to numerous practical problems; we will 
not consider any other numerical values for these two parameters. We now investigate another simple 

case, co * 0, ^ j • Then, we have 




3 Pr 


Substituting this expression into 9 t£?(^, T|), we obtain a sixteenth-order polynomial with respect to T|. 
This polynomial contains only even degrees, namely, 0, 2, 4, 6, 8, 10, 12, 14, and 16. It is possible to 
make sure (we always use Mathematica (ref. 21) to perform cumbersome transformations) that the coef- 
ficients of the above polynomial are positive for all co (co ^ 0) and for all Re; consequently, the corre- 
sponding sixteenth-order equation has no real roots. Therefore, the determinant (eq. (18)) has no other 
real zeros in this case as well. 


We have finally come to the most complicated case, which so far allows us only approximate inves- 
tigation. Let M% < ^ j . Then, we have to clarify whether the functions 9i CK^g^q ), q ) 


(j\ 

and/or 910(^3 (q), q) turn into zero forq within the range given in inequality (22). Both functions are 

actually of a general algebraic type (they contain non-integer powers), which means we have only a 
remote possibility of accurately (analytically) showing that they have no real roots, particularly because 
these functions depend on many parameters. At least at this point we are unable to construct the corre- 
sponding rigorous proof; therefore, we use the following graphical approach. 
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To start, we select some representative discrete set of the parameters involved. The range for the 
Mach number is known, so we simply choose a few points within this range. As for the Reynolds 
number, the representative values for the graphical tests we are conducting may be chosen to be about a 
few thousand. Indeed, we are not studying Stokes’ flows that correspond to very low Re. As for typical 
laminar solutions for the flows around an airfoil, they apparently cease to exist starting with Reynolds 
numbers at around a few thousand. Moreover, for turbulent flows with true molecular Reynolds num- 
bers of a few million, one can successfully model turbulence in the far field by introducing a new effec- 
tive value of the Reynolds number, which also appears to be around a few thousand. (See ref. 12.) 
Finally, recall that the periodicity of flow in time is caused by some external influence, and reference 7 
reports that the maximum effect (i.e., response) of such an influence corresponds to nondimensional fre- 
quencies of about 1 . Therefore, we will not consider frequencies much less than unity or much greater 
than unity. The upper bound for the band of frequencies originates from the numerics since we are 
going to pass from series (6) to the finite Fourier series while actually solving the problem on the com- 
puter. (See section 3.) 

We also note that the limits for r\ (see inequality (22)) do not depend on 
the sign of (0. Moreover, since £^(r|, M 0 , Re, Pr, y, |co|) = Re, Pr, y, -|oo|) , 

£g 2) (xi, M 0 , Re , Pr, y, |co|) = -{^(rj, M 0 , Re, Pr, y, — |co[ ) (eq. (20)), and all powers of £, and (co + ^) in 

T|) are even, it suffices to investigate the behavior of only one of the above functions for both 
positive and negative values of co. We do this by plotting the corresponding graphs for the following 
specific values of the parameters involved: co = ±0.5, ±1, ±10, and ±50; A/ 0 = 0.4 and 0.7; Re - 1000, 
2000, and 5000; and y and Pr are still 1 .4 and 0.72, respectively. The graphs drawn with the help of 
Mathematica (ref. 21) in different scales show that neither of the above curves intersects the real axis. 
(We do not present these plots here because they are not of interest except to show that the correspond- 
ing curve has no zeros). Relying on this approximate graphical investigation, we may expect that at least 
within some range of the parameters involved the symbol Q (eq. (16)) has no other real singular points, 
except for those that have already been found. 


We use an analogous graphical approach for the case M ft > ^ + p") *^ e ^ ave no P re " 


scribed range for r| in this case. However, it is clear that the asymptotics of the functions 
^2(4^’ 2 H r l)» T l) f° r large T| is T] 8 , so it suffices to study the behavior of the above functions only on 

some finite interval of rj. We used Mathematica (ref. 21) to plot the corresponding graphs for the same 
values of co. Re, y, and Pr as mentioned before and for A/ 0 = 0.8 and 0.9. The graphs (drawn in different 

scales for different r|-intervals, up to -10 5 < r\ < 10 5 ) show that neither of the curves has real zeros in 
this case as well. 


Summarizing, we conclude that at least for a certain reasonable range of the parameters involved, 
M 0 , Re, Pr, y, and co, we have justified the following proposition. 


Proposition 1: The symbol Q(^, ri) (eq. (16)) has only three real singular points on the (^, r\)-plane: 


(~co, 0), 


( coAf 


1 -M 


,0 

o ; 


, and 


co M 


1 + M 


,0 


0 


. For co = 0, these three points merge into one. 


To determine whether the inverse symbol Q _1 (^, T|) belongs to Lj oc (R 2 ), it suffices to investigate 
the behavior (integrability) of this matrix function near the three singularities. This investigation actu- 
ally means that we have to check integrability of each of the 16 elements of Q _1 (^> r\). These elements 
are given by (Q -1 ) 7 , / = 5; ■/(?, 1 < i , j < 4, where j are the corresponding cofactors. 
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Figure 2. Powers involved in the denominator QQ (black circles) and Newton’s diagram (dashed line) for QQ ; 

(^ = -to, i) = 0 ), co * 0 . 

Let us first concentrate on the singularity (£, = -to, r| = 0) for to * 0. We replace the above expres- 
sions for the elements of inverse symbol by their equivalents, (Q 1 ) y , = (8 ( . jQ)/(QQ) ( Q means 
complex conjugate), to make the denominator purely real. Since both the denominator QQ and the 
numerator 8 ( jQ are the sums of monomials of the type const ■ (to + £,)*T| / £, m (here const depends on M 0 , 
Re, Pr, y, and to, and k, l, m are nonnegative integers), then it would be sufficient to make sure that any 
expression of the sort 


|con.rr.(to + 5)V^ OT | 

QQ 

that originates from (8 ( jQ)/(QQ), 1 < i,j < 4, is integrable near (£ = -to, r| = 0). Since to * 0, then the 
factors do not contribute to the asymptotic behavior of expression (24) near (£, = —to, T] = 0), which is 
an essential difference in comparison with the case to = 0. (See the following discussion.) Therefore, we 
may investigate this asymptotic behavior by constructing Newton’s diagram (see ref. 22 by Walker) 
with respect to only two variables, to + £, and Tj. Namely, we show in figure 2 the set of points (k, l) that 
correspond to all monomials const ■ (to + %) k T\% m involved in QQ . The Newton diagram (ref. 22) is a 
lower part of the convex hull of the above set. The diagram is shown by the dashed line in figure 2. 
Those points ( k , /) which belong to the Newton diagram determine the asymptotic behavior of QQ near 

(5 = -a*. Tl = 0). 

More precisely, the asymptotic behavior of QQ near the singularity is determined not only by the 
lowest degree monomials (see Newton’s diagram in fig. 2) but may also depend on some higher order 
terms if the form 
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(25) 




16 

9 


Re 4 Pr 2 M$ 


q 8 + 


Re 2 M$ 


II + ?;T-5l;] ( " +5,v 


(which corresponds just to the previously mentioned lowest degree terms that constitute the Newton 
diagram) degenerates under some conditions. However, in this specific case the form (eq. (25)) is 


positive definite because 


4 1 \ 2 

3 + Pr) 


8 

3 Pr 


is positive for any Pr. Therefore, after some natural change 


of variables (see the following text) the asymptotic behavior of the denominator QQ becomes uniform 

with respect to the polar angle, which implies that while investigating the integrability of Q _1 (^, T|) one 
may simply neglect all the higher-order terms (black circles above the dashed line on fig. 2) and con- 
sider the expression 


|corm.(co + £)*V£ m | 

W’ 

instead of equation (24). Furthermore, we may only increase the ratio (eq. (26)) by neglecting the third 
term (~(co + £,) 2 r| 4 ) in equation (25). Indeed, it is easy to see that in doing so we only decrease the 
denominator but still preserve its positive definiteness. Finally, eliminate the factors for simplicity. 
We have already mentioned that £, m do not contribute to asymptotics near (£, = -co, T| = 0), (0 * 0. There- 
fore, to estimate the integrals, we may replace these factors by appropriate constants, e.g., 

\const • (to 4 - ^)^V^ m | < \const * ((0 + ^) T[ ll^l max 

(^/M £)(co + ^) 4 + ( 16/9) &/Re*Pr 2 M*ytf ' + £) 4 + ( 16/9) (l; 4 m / Re 4 Pr 2 M*)j]* 


where minimum (min) and maximum (max) are found on a sufficiently small neighborhood of 

(5 = r| = 0). 

Thus, we have reduced the original question of integrability of (8- jQ)/(QQ) to checking the inte- 
grability of the following function: 


\const • (co + ^,) k r\ l \ 
a(co + ^) 4 + br\% 


(a > 0; 6 > 0; k, l are nonnegative integers) 


(27) 


on some neighborhood of (£, = -co, r\ = 0). Because of the symmetry, it suffices to integrate func- 
tion (27) only on one quadrant, for example, co + £ > 0 and rj > 0. Moreover, since we are studying local 
integrability, we also introduce some upper limits for co + ^ and for T|, e.g., co + 2; < 1 and r| < 1 . Let us 
now change the variables Ja { co + £,) 2 = £ and Jbr | 4 = % and then proceed to the following integral: 


const 


l 1 

ii 

00 


£(*- l)/ 2 ^(/- 3)/4 

C 2 + x 2 


d^dx 


(28) 


Further, make another change of variables, from Cartesian (£, %) to polar (p, 0) coordinates, and for 
simplicity, truncate our rectangular domain, {0<£<l,0<x^ 1} -^ { C 2 + X 2 - 1 } » which obviously 
does not influence the result (integrable or not integrable). Finally, we obtain, instead of integral (28), 


[ + (*- l )/2 + (/ - 3)/4 


n/2 


const 


J (cos 0 )^ 1 ^ / 2 (sin 0 )^ 3 ) / 4 <i 0 <ip 


(29) 
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Figure 3. Black circles represent powers of the monomials in cofactors for = -co, T| = 0), co * 0. Gray area corresponds to 
integrability conditions (30). 


From formula (29) one can easily derive the conditions sufficient for the integral to exist. Namely, they 
are 


k - 1 
2 


+ 


1-3 

4 


>£ 


(30a) 


k- 1 

— > e - 1 (30b) 

1-3 

— >e-l (30c) 

where e is an arbitrarily small positive number. 

We now have to make sure that all conditions (30) are satisfied for all cofactors 8, ,, 1 < i,j < 4. 
First, we note that since k and / are always nonnegative integers, then two conditions (eqs. (30b) and 
(30c)) are met automatically. Then, to check the fulfillment of the third condition (eq. (30a)) one has to 
accurately calculate all monomials involved in all cofactors 8 (y , 1 < i,j < 4 and analyze the powers (k, l) 
for (to + fy k r\ l . This step was done with the help of Mathematica. (See ref. 21.) In figure 3, we have col- 
lected all the relevant powers (k, l ) for all cofactors 8 ( j, 1 < i,j < 4. We also show in figure 3 the range 
of powers (k, /) which satisfies conditions (30) (gray area). Using figure 3, one can easily conclude that 
all monomials involved satisfy all conditions (30). Therefore, the inverse symbol Q~ (^, r)) is abso- 
lutely integrable near the singular point (^ = -to, T| = 0) for 0) * 0. 
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Figure 4. Powers involved in the denominator QQ (black circles) and Newton’s diagram (dashed line) for QQ ; (£ - 0, r\ - 0), 
G>* 0. 


The integrability of Q _1 (£, r|) for co = 0 near the singular point (£ = 0, rj = 0) is investigated by the 
same method. We only note that since £ and co + ^ are now the same, both of them do contribute to the 
asymptotic behavior of Q _1 (^, r|) near (| = 0, r\ = 0). Therefore, the sets of monomials involved, as well 
as the Newton diagram, for QQ will differ noticeably from those relevant to the case co 0. Indeed, the 
asymptotic behavior of the denominator QQ near = 0, r\ = 0) is now determined by the following 
form (compare with expression (25)): 


,,<»» = (,+ ' J-Vs 

Q V Mft Mir 


,12 


M$Pr 2 Re 4 


_? ?_V6 T1 2 + + ! (\ + J_? _ 2 l 

M 4 Mir ^ M 4 Re 2 l Pr) Pr 


(31) 


which corresponds to the Newton diagram presented in figure 4. 

As in the case co * 0, the form A ( q (eq. (31)) also appears to be positive definite since all five coef- 
ficients in expression (31) are positive for all Re, Pr, and Mq < 1 . However, the Newton diagram shown 
in figure 4 consists of two straight intervals, whereas the one in figure 2 contains only one interval. This 
difference is essential because now each of the aforementioned two intervals (see the two-component 
dashed polygonal line in fig. 4) will determine its own domain of integrability for the expressions 


\const • ^ k T\ l \ 

4 (0) 

A Q 


(32) 
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on the ( k , /)-plane. Here k and / are the powers in the numerator of expression (32). Since the form A q 1 

is not simply positive definite, but all powers involved are even, and each coefficient in formula (31) is 
positive, we can find the corresponding domain of integrability on the (k, /)-plane independently for 
each of the two parts of the Newton diagram. (See fig. 4.) To do this for either part of the diagram, 
neglect those terms in the denominator which correspond to another part (in so doing, the denominator 
may only decrease). Then, formally divide both the numerator and the denominator by the common fac- 

tor q . Using the changes of variables analogous to those previously implemented, we come to the fol- 
lowing set of conditions sufficient for the integrability of function (32) near (£ = 0, r| = 0): 


1 

2 


1-7 
+ — — 


>e 


(33a) 


k- 1 
2 


>£- 1 


1-7 


>e- 1 


and 


k - 5 
2 


+ 


/- 1 
2 


>£ 


(33b) 

(33c) 


(34a) 


k - 5 
2 


> 8-1 


(34b) 


/- 1 
2 


> 8-1 


(34c) 


Note that three conditions (eqs. (33)) correspond to the upper part of the Newton diagram and three con- 
ditions (eqs. (34)) correspond to its lower part. (See fig. 4.) 

In figure 5, we show (with black circles) all powers (k, /) involved in all cofactors 5, p 1 <i,j< 4 
for the case to = 0. Gray areas on this figure correspond to the range of those coefficients (k, 1) which 
satisfy integrability conditions (33) and (34). Note that conditions (33c) and (34b) impose some addi- 
tional restrictions on / and k for the upper and lower components respectively of the Newton diagram in 
figure 5. We did not have such restrictions in the case to * 0. (See inequalities (30).) One can easily see 
from figure 5 that all elements of Q *, (5 J jQ)/{QQ), 1 <i,j< 4, are absolutely integrable near 
(^ = 0, t| = 0) in the case to = 0 as well. 


Finally, we only have to show that Q ! (^, q) is absolutely integrable on some neighborhood of each 


of the singular points 


f C0A/ 0 ^ 

,0 


1 - M 


0 


and 


J 




0 


1 +M. 


,0 


0 


for to * 0. If we simply ensure that Q , (£.'n)l is 


integrable on the same neighborhood, then the integrability of Q '(^. q) follows. To do this, first note 
that grad Q(E„ q) * 0 at either of these two points. Indeed, it is quite easy to see from equation (18) that 


^0 at both 


( to M 


1 -M 


, 0 and 


( coA/ 0 ^ 

,0 


o 


1 + M 


0 ) 


for all Mq< 1 . Then, refer to reference 23, wherein 


Vainberg proves exactly the same statement we need, namely, the integrability of |g -1 (i;, q)| on some 
neighborhood of an isolated real zero of the polynomial Q(q, q) when grad Q(lq, q) * 0 at this point. 
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0 1 k 


Figure 5. Black circles represent powers of the monomials in cofactors for = 0, T| = 0), co * 0. Light-gray area corresponds 
to integrability conditions (33). Middle-gray area corresponds to integrability conditions (34). Dark-gray area is common to 
both conditions (33) and (34). 

Thus, we can finally formulate the following proposition. 

Proposition 2: The inverse symbol T|) (eq< (16)) is absolutely integrable on any finite domain of 

R 2 , i.e., Q -1 (£,Tl)e Lj oc (R 2 ). 

In accordance with reference 20, proposition 2 immediately implies proposition 3. 

Proposition 3 (existence): The system (15) is solvable in Q' for any compactly supported f e L l (R 2 ); 

f in equation (15) is a Fourier transform off (eq. (14b)). 

The solution to the AP (eqs. (12) and (13)) that we are looking for may generally be found by 
means of the inverse Fourier transform (again, the superscript n is omitted below), 


u(x, y) = L J J ri)e‘^ + '^4rfri (35) 

— oo — oo 

* * V 1 a v 

Using the brief notation, we may rewrite formula (35) as u = (u) = (Q~*f) . However, in doing so 

we still do not know whether the function u(x, y) (eq. (35)) satisfies boundary condition (13). Let us 
first prove the following proposition. 

Proposition 4 (uniqueness): If the solution u of system (12) satisfies the boundary condition (13), then 
it is unique in the class of distributions vanishing at infinity. 
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Indeed, any function u that solves system (12) is actually an inverse Fourier transform of some solution 

to system (15), u = (u) v . In turn, any distribution ue g' that solves system (15) (see formula (17)) 
should coincide with the regular function Q (q. r| )f( q, r|) everywhere on R 2 except at the three singu- 
lar points of Q(^, T|) (since f(£, r|) has no singular points). Therefore, any other solution to system (15) 
may differ from u only by a distribution with the support belonging to the three-point set 


(-to, 0), 


( co M 


1 -M 


,0 

0 J 


coM 


1 + M 


,0 

0 ) 


■ Such a distribution may be only a finite sum of 8-functions and 


their derivatives. (See ref. 20.) Therefore, if u = (u) v vanishes at infinity, then any other solution to 

system (12) will differ from u by an inverse Fourier transform of a finite sum of 8-functions and their 
derivatives, and, consequently, it will not vanish at infinity since Fourier transforms of 8-functions and 
their derivatives are polynomials. (See ref. 20.) Thus, proposition 4 is justified. 


Let us now select a finite ball B where 


B c R 2 , B - ED (-CD, 0), 


c oM 


1 -M 


,0 


0 


(0 M 0 
1 + M { 


,0 

0 / 


(e>0) 


and construct a partition of unity, 1 = g 5 + g^, where the functions g B and g- are infinitely smooth 

and bounded on R . The function g B is identically zero outside the ball B + e, therefore, g- is identi- 
cally zero inside the ball B - e. Note that such functions always exist. (See, e.g., ref. 20.) Obviously, 

_.i v iv ; v 

u = (Q *f) = (Q ‘ggf) +(Q 'g-f) • We will separately analyze each term on the right-hand side 
of the above sum. First, it is clear that Q“'g B fe L l (R 2 ) because f is bounded and Q 1 e L} oc (R 2 ). 

* V , £ V 

Therefore, (Q _1 g fi f) — >0 while Jx 2 + y 2 — »oo. For the second term (Q -1 g-f) we cannot yet con- 

struct a general proof of its decay at infinity. The difficulties here arise from the fact that 
Q -1 g L} oc (R 2 ) but Q -1 € L l (R 2 ), i.e., it is not absolutely integrable near infinity. Therefore, a general 

proof may require an appropriate regularization of the corresponding oscillatory integral. However, we 
retain this question for a future investigation. For the time being, we can formulate the following two 
statements. Each will address the vanishing of the solution at infinity for some particular case (or in a 
weaker formulation). 


First, assume that f e L 2 (R 2 ), which is actually not restrictive for our purposes. Then, f e L 2 (R 2 ) 
(we may treat the Fourier transform here in the sense of Plancherel, ref. 24). As mentioned before, 

Q -1 e L l (R 2 )-, however, Q“’g„ can be shown to be bounded on R 2 . Therefore, Q -1 g-fe L 2 (R 2 ) , 

tf B 

: v 

which immediately yields (Q“'g-f) e L 2 (R 2 ) . Thus, in this case the solution u to system (12) is rep- 
resented as a sum of two terms, u (1) + u (2) , where u (l) -> 0 while Jx 2 + y 2 -> °° (true vanishing in the 
sense of boundary conditions (13) and u (2) e L 2 (R 2 ), which may be treated as a “generalized decay”. 
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We also note here that the statement on uniqueness proven in proposition 4 also applies to the functions 
from L 2 (R 2 ) since the polynomials obviously do not belong to L 2 (R 2 ). 

Second, if we impose some additional restrictions on f , namely, if we require that f be sufficiently 

smooth on R 2 so that fe L l (R 2 ), then we obtain a true decay for the second term as well, 

(Q-'gJ) V ^0 while Jx 2 + y 2 — » . Therefore, for a more particular class of the right-hand sides we 

B 

may affirm that the problem (eqs. (12) and (13)) is uniquely solvable in g\ We note that for some other 
cases (see ref. 9) such a restriction of the class of admissible right-hand sides does not influence the con- 
struction of a DPM-based numerical algorithm. We will not rigorously formulate and prove this state- 
ment for the specific case currently under study. However, we expect that this property does take place. 
These expectations are based on the numerical experience. (See refs. 8, 1 1, and 12.) 


2.3. Truncation of Linearized Problem 

As mentioned in section 2.1, we are not going to directly solve the problem (eqs. (8) — (10)). Instead, 
we will implement some additional truncation and further solve only a new finite substitute for the lin- 
earized problem. Since the problem (eqs. (8)— ( 1 0)) will be solved by means of the DPM, we must con- 
struct an equivalent finite substitute for the auxiliary problem (eqs. (12) and (13)). Moreover, the same 
finite substitute for the AP (eqs. (12) and (13)) will be used for calculating the operator , which pro- 
vides boundary data for the problem (eqs. (8)-(10)). (See section 2.1.) In this section, we construct the 
finite substitute for the AP introducing some additional assumptions in regard to both the smoothness of 
the solution we are looking for as well as the rate of its decrease at infinity. This is done in order to sim- 
plify the presentation and to avoid unnecessary complications that are not essential for the purpose of 
constructing the numerical algorithm. We hope to provide a more rigorous analysis of the approach 
described here in a forthcoming paper. 

For reasons of numerical convenience and effectiveness, we will use a different method for calcu- 
lating the solution of the AP, rather than the one from section 2.2. Using this new solution technique, we 
will equivalently reformulate the AP on a new finite domain. Namely, let us again take the Fourier 
transform of both sides of system (12); however, now we do so only in one Cartesian direction, y (com- 
pare with eqs. (14)), 


u(*,ti) = ~^= f u(x, y)e~ ir *ydy (36a) 

j2n J 

— oo 
oo 

f(jr, ri) = — = f f (x, y)e~'^> dy (36b) 

J2n J 

(Again, we drop the subscript n hereafter in this section to simplify the notation. Moreover, we retain 

here the same notation, u and f , as in section 2.2; however, the left-hand sides of expressions (14) and 
(36) are obviously not the same.) Then, we obtain the following family of systems of ordinary differen- 
tial equations (ODE’s): 


d u(x, T| ) 
dx 


+ Q(ri)u(x, r|) = f(x, T|) 


(37) 
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where 


and 


QOi) = D 


- n-l 


(CO + 


Re 


*T 1 
0 


(CO + 


4 T) 2 
3 R~e 


0 

0 


1 


ICO 

0 

0 



0 


. co + jni _uo IlL_ 

RePrMl 


f(-r> T|) = D ^(jc, rj ) 


(the matrix D is defined in formula (36)). The family (eq. (37)) is parameterized by the continuous vari- 
able r|, < r| < oo, and * is an independent variable. Recall that the solution u(x, y) we are going to 

calculate should vanish at infinity. (See boundary conditions (13).) Consequently, we will generally 
impose the following boundary condition on the solution of system (37): 


ujc, q -> 0 as |jc| -» °° (38a) 

However, in particular cases (see the following discussion and ref. 8 for details) the condition 
(eq. (38a)) may appear too restrictive — namely, the cases when Q(q) has purely imaginary (or zero) 
eigenvalues. Therefore, for some selected values of co and r| we will only require 

|ujc, T|| < const as |jr|-»°o (38b) 

Note that we do not consider solutions that grow polynomially, the latter solutions correspond to the 
case when Q(q) has multiple purely imaginary eigenvalues and does not have a basis composed of 
eigenvectors. 

Once we are able to find (for every q) a solution to system (37) that would satisfy boundary condi- 
tion (38) at infinity, then the solution to the AP (eqs. (12) and (13)) can be restored by means of a one- 
dimensional inverse Fourier transform, 


u(*,y) = 


4 = f u(jc,q)e' t l>rf n 

V2tc j 


(39) 


Let us designate the inverse operator for the one-dimensional problem (eqs. (37) and (38)) by G^.(q). 
That is, the solution u(x, q) to this problem is given by 

u(*,q) = G x (q)f(*,q) (40) 

The operator G ( (q) is obviously linear. Combining formulas (36), (39), and (40), we obtain the follow- 
ing formula for the solution of the AP (eqs. (12) and (13)): 


u(jc,y) = ^ J G v (q)| f(jr, (41a) 
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Now, we will show how one can pass from the AP (eqs. (12) and (13)) to the new AP formulated on 

Y < < Y 
2 y 2 


Y y 1 

the strip {-oo<^<oo}x<j — - < y < - > and periodic in the y direction, with Y being the value of the 


period. In doing so, we expect that when the period Y grows, Y -4 the solution to the new periodic 
AP will uniformly converge to the solution of the original AP (eqs. (12) and (13)) on any strip 
{ — oo <jc<oo}x { — y < y < y } where y is fixed and always less than Y/2. We note that the same approach 
was used in reference 8 for the steady-state problems. 


Hereafter, we assume that all functions involved are defined on the infinite strip 
f y Y } 

{ _oo < < oo } x < -- < y< — >. The width of the strip Y is initially supposed to be greater than the 

diameter of suppf . (Later we will consider the limit Y — > oo.) We assume periodicity of the solution to 
the new AP in the y direction. Then the solution that vanishes as \x\ — > °° is given by 


uy(*, y) 




k — — « 


Y/2 


*{ Y JY 


J ?(•*> 


s)e 


.Ink, 


y) 


ds 


-Y/2 


(41b) 


In formula (41b), we use the Fourier series of a periodic function instead of the Fourier integral used in 
formula (41a). Our aim is to estimate |u(jc, y) - uy(x, y)| from above on a finite (fixed) interval (-y, y), 
y < y/2 , uniformly with respect to x . Let us introduce a uniform mesh in T|, where h ^ = 2k!Y is the 

mesh size, and designate r|* = kh^ k = 0, ±1, ±2 Let us then fix some interval (-JZ, Jf); we will 

always choose h ^ (and consequently Y) so that + 1/2), A" being an integer. Then, 


u(*,y)-uy(*, y) = — 


\_ 

2n 


Y/2 


Jg^ti) £ VW J f \x,s)e iX[ * (s - y) ds 


| — oo 

— oo 


k = ~oo 

<± 

2n 

II 

1 

+ — 

2n 


j...- x - 

I X 

-A k = -K 

|t|| >a W>K 


-Y/2 


-s'v&h. 


Let us separately estimate each of the two terms (the first one corresponds to the finite interval, and the 
second one corresponds to the complementary infinite interval). 


k-K 



k = -K 
k = K I 



k = -K 


(tf+1/2)/^ oo Y/2 

J G^ri) J i(x,s)e- i <*-y'> 1 ldsdr\-h 1i G x (r\ k ) J f \x,s)e~ i ^ k(s ~ y) ds 

{K-\/2)h rx -oo -Y/2 

'( K+l/2)h ,, 

J G x (ti) | {(x,s)e-‘^-y^dsdn-h J) G x m k ) J f(x,s)e^ k(s ~ y) ds 

(tf-1/2 )* 


+ /? 


n 


G^T)*) J f(*, s)e ,T1 * ( ^ y) ds 

|j| > Y/2 
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Clearly, the right-hand side of this inequality is actually the sum of errors of the quadrature formula of 

a r 00 a 

rectangles for the function u(x,T\)e‘^>’ = G^(r|) f(jt, ( se e formula (40)) on elemen- 

* —00 

tary segments of the kind [(*- 1/2 )h n , (k+ 1/2)/^], k = -K, Indeed, for each k,k = -K , ..., K, the 
third term that corresponds to the integration over \s\ > Y/2 turns into zero for sufficiently large Fs 
since f(x, .v) is compactly supported. Therefore, one can obtain the following estimate: 


— | |, < const ■ h?:A max 
27C 1 n xeR 

T1 e (-j?, A) 


a 2 


< const ■ max 

n xe R 

p 6 (-JI A) LI 


U^T])^ 

+ 2|y| 


3rj 2 

a 2 u(x, q) 


dri 2 


du(*. T|) 

dii 


+ y 


! |u(x,n)| 


= (c, +c 2 \y\ +C i y 2 )h^ (c,,c 2 ,c 3 > 0) 


Note that if we initially assume that the solution u(x, y ) decreases at infinity sufficiently fast, then the 

differentiability of its Fourier transform u(x, n) (see the right-hand side of the previous inequality) fol- 
lows directly. 

For the second expression, we obtain 


— h 2 < — 

2n 12 2 tc 



OO 


Y/2 


! 

G v (r|) J f(jt, 


G^) f f {x,s)e~ i ^ s ~y)ds 

J 


hi >* 

—00 

1*1 >* 

-Y/2 



Y/2 


Let us replace the integration limits J in the second term on the right-hand side of this inequality by 
~ -Y/2 

J , as was done when estimating |*| j. Then, 


±|.| <_L 

2 71 ' 2_ 27t 


J |u(x, r|)|dr| + ft tl |u(x,Ti t )| 
hl>^ W >* 


Additionally we assume that the solution we are looking for has two absolutely integrable derivatives. 
Then its Fourier transform decreases faster than |t)|~ 2 and the previous inequality straightly implies 

< c ->o> 

Combining the two obtained estimates, one easily gets 


|u(x, y) - u Y (x, y) \ < c Q /i 2 A + 


where c 0 ^jnax ^ (c, + c 2 |y| + c 3 y 2 ), c 0 > 0. Clearly, all constants involved in the foregoing esti- 
mates depend, generally speaking, on the specific nonperiodic function u(jc, y) that we approximate by 
the periodic functions u K (x, y). 
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Now let e be an arbitrary positive number. We will choose sufficiently large Y e (i.e., sufficiently 

small h ) so that the following inequality 

*»£ 


h 2 A + — <e 




A 


(42) 


is satisfied for all Y > Y . In other words, we require that for a prescribed e and for any h n < /i„ ine- 

quality (42) has real positive solutions A of the special kind. A - {K + 1/2)/^ ( K being an integer). The 
latter requirement is always met if, e.g., the distance between the real roots of the quadratic equation 
c 0 h 2 A 2 -eA + c 4 = 0 is greater than h This, in turn, yields the inequality e 2 -4 c 0 c 4 h 2 >0 

for hr*. This inequality is obviously satisfied for any 0 < h < h , where h e R is a unique positive 
I M Me Me 

root of the equation e 2 - 4c 0 c 4 h 2 - = 0. Since the fulfillment of inequality (42) is sufficient for 

the estimate 


u(x, y) - u Y (x, y) <e 


(43) 


to be true, then we have shown that for all e > 0 one can always find a sufficiently large period Y z so 
that for any Y> Y E the absolute value of the discrepancy between the nonperiodic solution u(jc, y) 
and its periodic approximation uy(x, y), |u(jc, y) -Uy(jt, y)|, does not exceed £ for all x and for all 
-y <y<y. 


Thus, we have reduced the original AP (eqs. (12) and (13)) to the new AP formulated on the strip 


{_oo<jc<oo}x 



. In section 3, we show that we will only need to know the solution of the 


AP in some neighborhood of suppf , therefore, the approximation of the nonperiodic function u(jc, y) by 
a periodic one, u y(x, y ), only on a finite interval -y<y<y is sufficient for our purposes. Let us now 

f Y f] 

show how to pass from the domain {-o°<jc<«>}xj — <y<-k which is still infinite, to a truly finite 
domain for the new AP. 


Instead of 


{-oo < X < oo } x 


Y Y 
-- < y < - 
2 y 2 


let us now consider a rectangular domain 


D £ = (0, X) x (-T/2, Y/2). (See fig. 1.) This new domain D £ should completely contain T and Tj. 

We will reformulate the new AP so that its solution will be determined only on this finite domain D £ 
and will coincide there with the corresponding fragment of the solution found on 

f y y } 

{-oo <jt<°o}x^--<y<-| before the reformulation. As previously mentioned, we only need to cal- 
culate the solution to the AP in some neighborhood of D in . Therefore, we are always able to choose an 
appropriate X and Y so that this neighborhood belongs to D £ , and consequently we only need to con- 
struct special boundary conditions at the lines x = 0 and x = X so that the reformulated new AP being 

solved on D ^ is equivalent to the periodic AP on the strip {-°o<*<oo}x < y < ^ j described ear- 

lier in this section. These boundary conditions at ^ = 0 and x = X will be set separately for each 
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wavenumber k y k = 0, ±1 , ±2, . . (see formula (41b)) involved in the Fourier representation of the func- 
tion uy(x, y) . Namely, for each k, k = 0, ±1, ±2, ..., we require that the corresponding Fourier mode, 

u (*> q*) - u (*> 2nk/Y), meets boundary condition (38) at infinity. To exactly transfer condition (38) 

from infinity to the finite boundaries x = 0 and x = X, we use the following consideration. Since sys- 
tem (37) consists of ODE’s with constant coefficients and since it is homogeneous outside (0, X) (recall 

that suppf(*, y) c D in , and consequently suppf(jc, q) c (0, X) for all q), then it obviously has four lin- 
early independent eigensolutions (in the region of homogeneity). Depending on the structure of the set 
of eigenvalues of the matrix Q(q), these eigensolutions may either increase or decrease (eq. (38a)) 
exponentially, or they may oscillate (eq. (38b)) while x — > °° and while x — >— «*. As previously men- 
tioned, we do not consider the last possible case when Q(q) has multiple purely imaginary (or zero) 
eigenvalues and does not have a basis composed of eigenvectors, which leads to polynomially growing 
solutions. Sometimes one can analytically make sure that this case really does not take place. For exam- 
ple, we do this in section 3 in the discrete formulation for some particular values of co and r|. In other sit- 
uations, this question may require some additional numerical investigation as in reference 8. At any rate, 
to satisfy boundary condition (38), we must prohibit at x = 0 all solutions that do not decrease to the left 
(i.e., as * — » -oo ), and prohibit alx = X all solutions that increase to the right (i.e., as x —> °o ). The rea- 
son for this asymmetry was mentioned before: once we have purely imaginary (or zero) eigenvalues of 
Q(r|) and, consequently, oscillating or constant-in-space solutions (see formula (38b)), then we cannot 
always prohibit at both ends of the interval (0, X) all modes that do not decrease in the corresponding 

direction. However, it should not affect the result since the final solution we are looking for ( u(x , y)) 
decreases at infinity. (See subsection 2.2.) Moreover, we have proven in reference 8 that once we have a 
selected nondecreasing mode in Fourier representation of the solution, then after the inverse Fourier 
transform the entire solution will nevertheless decrease. Therefore, we can take into account selected 
nondecreasing modes (if any) by simply admitting them at one of the two boundaries, * = 0 or ^ = X. (In 
case we do not do this, the problem may appear overdetermined.) It seems more natural to admit the 
nondecreasing Fourier modes (if any) at the downstream boundary x = X. (See ref. 8.) 

Now, we calculate the eigenvalues X,{q*), r= 1, ..., 4, for the matrix Q(q*). Those eigensolutions 
that increase to the right correspond to eigenvalues 9tA. r <0, and those eigensolutions that do not 
decrease to the left correspond to eigenvalues > 0. Therefore, the following boundary conditions at 
x = 0 and x = X may be considered to provide an exact transfer of boundary condition (38) from infinity: 

S”(n*)»(0,Tl*) = 0 (k = 0, ±1, ±2, ...) (44a) 

S + (q*)u(X,q*) = 0 (A = 0, ±1, ±2, ... ) (44b) 

Here S~(q*) and S + (q*) are the special rank-deficient 4x4 matrices that depend on Q(q*), with their 
ranks equal to the numbers of eigenvalues ^(q^) with nonnegative and negative real parts, respectively. 
These matrices are given by 


s (v = n ( 45a ) 

9U r (Tl*)<0 


s+ oi*) = f[ (QOW-MW < 45b > 

SRMn*)^o 
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Here I is an identity matrix and products in formulas (45) are calculated in accordance with the multi- 
plicities of the eigenvalues. Analogous conditions will be used in section 3 while dealing with the finite- 
difference formulation of the AP. 

Thus, the formulation of the new finite AP is now complete. Namely, we have to solve equa- 
tions (12) for the compactly supported right-hand side f (suppf aD in ) on the domain D £ (see fig. 1) 
with the periodicity boundary conditions in the y direction (Y being the value of the period) and with 
boundary conditions (44a) and (45a) at x = 0 and (44b) and (45b) at x = X. In the next section, we pro- 
ceed to the finite-difference formulation of the problem and describe the numerical algorithm for setting 
the global DPM-based ABC’s. 

3. Numerical Method 

3.1. Finite-Difference Scheme 

Let us introduce a uniform Cartesian grid in x [0, 7], with k x , 
grid in x , y, and t directions, respectively. We designate this grid fA£° , 

*C° T = {(x m ,y j ,t l ) = {rnh x ,jh y -Y/2,lT)\h x ,h r 
m = 0, 1, .... M, M = X/h x \ j = 0, 1, .... 2 J+ 1, 2 J 
1 = 0,1, 2L+ 1,2 L+ 1 = T/x} 

We will construct a second-order finite-difference approximation 
grid 5\£° r (see formula (46)) using the stencil shown in figure 6. 


h y , and x being the sizes of the 
x> 0; 

+ 1 = Y/hy, 

(46) 

of the system (eq. (3a)) on the 



Figure 6. Stencil. 
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Namely, we use the first-order differences in the x and t directions and second-order central differ- 
ences in the y direction, and we center the scheme with respect to the point (m + 1/2,;, / + 1/2), which 
yields 


1 /u ^ + J — u ^ _ u ^ 

C| m ’ - Ufn, j + Ufn — / Uft 

2 \ x t 


">±WV-D 


/|| ^ + • — U^ + ? 1 . — 

m + 1 , j m, j m + 1 , j m, j 


) 


f n l 


./+ 1 


+ 4 F 


m, j + 1 u m, y - 1 + U m+l,y+l U m+l,y-l + U m[ j + 1 U wT j- 1 U m + l,y + 1 ~ u m + l,y- 1 


2h . 


2A. 


2A. 


,/+l ji/+l 

+ 


- u 

2h 


/ + l 






. + u 


./ 


- 2u 


m » 7 + | ^ 7 u m, j — I u m + 1,7 + 1 z ' u m + 1 , j + + 1 , j - 1 


V 


,/+ 1 

V 7 + 1 


2ui + ? + u /+ 1 


./ + l 


■1 1 t * _ 9 iW + 1 4. || J + 1 A 

~m, j tn, j - \ ^ /w + 1,;+1 m+1,7 U m+ 1,7-1 

ft , 2 , 


(47) 


The finite-difference scheme (47) is written for the nodes (m, j, /), m = 0, 1, ..., M- \J = 0, 1, 27, 

/ = 0, 1, . . 2 L with the assumption that we later impose periodicity boundary conditions in time as well 
as in the y direction. Note that the stability of the scheme of type (47) was examined for the model scalar 
equation 


du + du du _ 1 d 2 u 

dt dx dy R e 2 dy 2 


(48) 


It turns out that the corresponding finite-difference scheme for equation (48) is unconditionally stable in 
the von Neumann sense. 


Then, using the periodicity conditions (compare with formula (4)) 

U m,y = u mf7‘ (« = 0, 1 , .... M; j = 0, 1, .... 27 + 1 ) 
we implement a discrete Fourier transform in time (compare with formulas (7) and (6)), 


(49) 


“ 2 L+ 1 S U ‘m,j e 
l = 0 

L 

u Lj = X 


2 L 


-ini T 


2jc 


(m — 0, 1, M\ j = 0, 1, . 2 J + 1 ; n — — L, . L) (50) 


ini x 


2k 


(m = 0, 1, M ; 7 = 0, 1 , ...,27+ 1; / = 0, 1, .... 2L) 


(51) 


n=-L 


and instead of system (47) obtain 


, n _ n c 

+ u m + Ij) + C n D ^ - + ^-F 


r:.n 


u m, j + 1 u /n, j - 1 + u m + 1, j + 1 ~ u m + \,j- 1 


2 h 


2h. 


+ y+1 ^ u m, ) + u m, j- 1 u m +1,7+1 ~ 2u^ + 1, j + U^, + 1, j - 1 




*, 2 


= 0 


(52) 


Here s n = 2/sinQnty j/x ; c„ = cosQnlyl; n = -L,...,L; m = 0, 1, .... M - 1 ; and 
y = 0, 1, ...,27. 
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The finite-difference system (52) is a discrete analogue of the continuous system (8) on the two- 
dimensional grid ?\£ 0 , 

= Ux m ,y i ) = (mh x ,jh y -Y/2)\h x ,h y >0; 

m = 0, 1, = X/h x ; j = 0, 1, ...,27+1,27+1 = Y/h y } (53) 

In formula (53) h x and h y are the same as in formula (46). 

We also note that once t — > 0, then s n — > i2nn/T = /co n (see section 2) and — > 1 . 


3.2. Difference Auxiliary Problem 

Let us construct another Cartesian grid in D £ , 

5W° = { (* m + 1/2 , y t ) ■((».+ 1 /2)h x , jh y - JV2)|^, h y > 0; 
w = 0, 1, ..., Af - 1, M = AY^; j = 0, 1 27+1,27+1 = f//^} (54) 

The grid sizes h x and h y are the same as before. The difference AP is formulated for the inhomogeneous 
counterpart of system (52) with a certain compactly supported right-hand side. The unknowns for the 
difference AP are defined on the grid (see formula (53)), and the right-hand side is defined on the 
grid SW°. (See formula (54).) In doing so, we obviously have the second order of approximation. We 
will define the specific right-hand side for the AP, t„ + \/ 2 j , m = 0, M j = 0, 1, ..., 27, in 
section 3.3. As for now, assuming that this right-hand side is already known (suppf^+ 1 / 2 , j cD ( „),we 
provide an exact formulation of the difference AP and describe an effective algorithm for its numerical 
solution. 

In accordance with the results of section 2.3, we impose the periodicity boundary conditions in the 
y direction, 


4.0 = u m, 27 + 1 (m = 0, 1, M) 

<-l = <2 j (m = 0, 1, M) 

Then, we implement a discrete Fourier transform (compare with formulas (36)), 


u m, k 


2J + 1 


2J LL 2n 

~ lk J h y~yT 




(m = 0, 1, ..., M; k = -7, ..., 7) 


7 = 0 
27 


.... 2n 
-tkjh — 


fm + l/2Y = 2 tVtX^ + 1/ 2-^ (m = 0, 1, ..., A/ - 1 ; k = -7, ..., 7) 

7 = 0 

and obtain, instead of the inhomogeneous counterpart to system (52), 


(55a) 

(55b) 

(56a) 

(56b) 


^ + B”u m ^ — fm+\/2,k 

(m = 0, 1, M - \ ; k = -J, ...,7) 

(57) 

where the 4 x 4 matrices A£ and B£ are given by 



c„ 

cr. ct. 


A" = 2?C + -^D + 

k 2 h x 

-^F + -^H 
2 2 

(58a) 

V C„ 

CT, C Jl 


B *" = 2 C -r D + 

_^Jf + -^h 
2 2 

(58b) 
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Here r k - ^h y , t k = -4sin 2 Q&7i v ,^ j//i 2 , and C, D, F, and H are defined in for- 

mula (3b). For each wavenumber k, k = —J, ... J, system (57) is composed of ordinary difference equa- 
tions, and it is a discrete analogue of system (37). To find a solution to the difference AP, we will have 
to solve system (57) for all k, k--J, ... J. However, the formulation of the difference AP is still incom- 
plete. To complete it, we have to set some boundary conditions at m = 0 and m = M (as was done at 
x = 0 and x = X for the continuous case in section 2.3). These boundary conditions should guarantee the 
desirable far-field behavior of the solution (i.e., decay at infinity). They will be formulated separately 
for each wavenumber k,k = -J,...J, i.e., the system (eq. (57)) will be supplemented for each k by some 
boundary conditions at w = 0 and m = M. The idea for constructing these boundary conditions in the 
discrete case is analogous to the one implemented in constructing boundary conditions (44) and (45) for 
continuous system (37). Namely, when formally considered on an infinite one-dimensional mesh, 

< m < °°, system (57) obviously becomes homogeneous at least for m > M and m < 0. The homoge- 
neous system has four linearly independent eigensolutions: those that correspond to ||i"(£)| < 1 

decrease to the right (i.e., as m — > °°); those that correspond to |p"(&)| > 1 decrease to the left (i.e., as 
m —> -°°); and those that correspond to |p"(k)| = 1 have either constant or oscillatory behavior. Here, 

H?(k), r= 1, ... ,4, are the eigenvalues of the matrix Q£ d =^(A£) -1 B£. Let us note that while calculat- 
ing the eigenvalues (i r ( A j for the stationary Navier-Stokes equations (ref. 8) (the eigenvalues are calcu- 
lated numerically using standard NAG subroutines), we have found that for all specific sets of the 
parameters involved (i.e., grid sizes h x and h y and hydrodynamic parameters A/ 0 , Re, Pr, and y) the abso- 
lute values of eigenvalues were never equal to unity except for the case of zero wavenumber, k = 0. For 
k = 0, we have obtained a multiple eigenvalue |p(0)| = 1 . (See ref. 8.) However, even in this case the 
system matrix still has a basis composed of eigenvectors, which provides us with the reason for not con- 
sidering the polynomially growing solutions in reference 8. For system (57), we also have a particular 
case when the eigenvalues of the system matrix become equal to unity in absolute value. Namely, it is 

easy to see from formula (58) that Qg = (A$) _1 B$ = -I (identity matrix). Obviously, Q$ has four lin- 
early independent eigenvectors; therefore, we do not have polynomially growing solutions in this case 
as well. As for other values of k and n, a numerical check (as was done in ref. 8) will always be neces- 
sary to determine whether the eigenvalues ||i"(£)| = 1 exist. If such eigenvalues do exist, a check is 

also necessary to determine what their multiplicities are and if there is a basis composed of eigen- 
vectors. Relying on our previous experience (ref. 8), we assume that while solving system (57), we can 

restrict ourselves by considering only these two cases: |p?(k)| * 1 and |p"(k)| = 1 with the full system 
of eigenvectors. Nontrivial Jordan blocks (of order more than 1) for |p"(k)| = 1 are excluded from con- 
sideration. Note, if the basis composed of eigenvectors does exist for |p"(it)| = 1 , then system (57) will 

be treated exactly in the same way as in the case |p"(k)| * 1 (the only difference is that the stability con- 
stant becomes proportional to M). 

Returning to the question of setting the boundary conditions for system (57) at m = 0 and m = M,we 
require that, analogous to the continuous case (see section 2.3), boundary conditions at m = 0 should 
prohibit all modes that do not decrease to the left (i.e., asm-) -°o) and boundary conditions at m = M 
should prohibit all modes that increase to the right (i.e., asm^ °°). Therefore, we may represent the 
desirable boundary conditions in the form of matrix relations (compare with formulas (44)), 
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s" (*)£&* = 0 


(k = j) 


(59a) 


S* + (*)uJf,t = 0 (k = -J,...,J) (59b) 

where 

S"“(t) = f] (QJ-|i?(*)I) (60a) 

|n?(*)| > 1 

s n+ (k) = fT (Q n k -\i?{k) I) (60b) 

K”(*)| S 1 

(compare with formulas (45)). 

Thus, the final formulation of the difference AP is the following. One should solve the inhomoge- 
neous counterpart to system (52) in D £ on the grid (see formula (53)), where the right-hand side 
f£ + 1/2, j (suppf„ + 1/2,7 C ^in) ls s P ec ifi e( 3 on the grid (see formula (54)), with periodicity 
boundary conditions (55) in the y direction and boundary conditions (59a) and (60a) at the line m - 0 
and (59b) and (60b) at the line m-M. 

To solve the difference AP, we implement the following numerical procedure. First, apply discrete 
Fourier transform (56) to both sides of the finite-difference system, then solve the system of ordinary 
difference equations (57) with boundary conditions (59) for each wavenumber k, k = -7, .../, and 
finally restore the solution by means of the inverse Fourier transform, 

k - J A ik . h 2k 

U n m ,j = ^u^ k e JyY (m = 0, 1, j = 0, 1, 2J) (61) 

k = -J 

The type of boundary conditions (59) (which are imposed separately for each wavenumber k) makes the 
choice of numerical method most relevant. An effective algorithm for solving one-dimensional prob- 
lems (eqs. (57) and (59)) is delineated in our work (ref. 25). We do not reproduce the corresponding 
results here, we only note that this algorithm may be thought of as a version of the well-known succes- 
sive substitution technique but without its “inverse” or “resolving” part. The computational cost of the 
numerical procedure in reference 25 as applied to solving the problem (eqs. (57) and (59)) is 0{M) 
operations (for each k, k = — .../). 

Let us now briefly describe the concept of convergence for the solutions of the difference AP. 
According to section 2.3, we approximate the nonperiodic solution by a periodic one on a finite interval 
-y < y < y when the period Y grows, T — » » . In its own turn, an approximate solution to the periodic 
problem is found by a finite-difference method on the grid with sizes h x and h y . Therefore, we will con- 
sider (uniform) convergence of the periodic difference solution (i.e., solution of the difference AP) to 
the nonperiodic continuous solution (i.e., to the solution of the original continuous AP) only on a finite 
rectangle (0, X) x (-y, y ) (this rectangle should be large enough to contain at least Tj) rather than on the 
whole domain of the difference AP. Moreover, we will consider this convergence not only when the 
grid size vanishes but also when the period Y synchronously increases, i.e., as ( h x , h y , Y) — > (0, 0, «>) . 
Of course, the rate of decrease for the grid sizes h x and h y and the rate of increase for the period Y are 
not independent; some estimates connecting these rates can be found in reference 8. Furthermore, some 
numerical experiments from reference 8 show that the presented construction of the difference AP does 
ensure the convergence of its solution to the solution of the continuous AP in the sense just described. 
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3.3. Computation of ABC’s 


In accordance with section 2.1, to set the ABC’s we need to know the following data: 


Uy, 




Here, £ is the normal to T. When integrating the Navier-Stokes equations step-by-step in time, we 


assume that 


f 

U V , ^ 


is provided from inside D in \ then we use these data to restore uj , which enables 


us to advance the next time step. However, as we carry out our analysis in the Fourier space, we cannot 


consider 


Uy, 


dC 


u" ^ 

Uv ’ ac 


as the actual values obtained inside the computational domain. To get 

. Without loss of generality, we 


we first have to calculate the Fourier transform of the function 
may always think that the latter is specified at the following points: 


( 

“ v ’ 


v x { t/| / = 0, 2L+ 1, t(2L + 1 ) = T} (62) 

Of course, actual discretization in time for the Navier-Stokes equations inside Dj n should not necessar- 
ily coincide with the one used for the solution of the exterior linearized problem. (See formula (46)). 
However, we may always use some interpolation in time to obtain the boundary data on a uniform mesh 
with respect to t (eq. (62)), which is convenient for further consideration. Hereafter, we simply assume 
that this interpolation (which is one-dimensional in time and of sufficiently high order) has already been 
implemented for each node v, if necessary. 


Another important issue related to the step-by-step integration in time is that the function 


u v- 


which provides the boundary data, is not necessarily time-periodic until we achieve a true oscillatory 
regime. However, for the purpose of constructing the ABC’s, we will propose some generalized treat- 
ment of the boundary data as being already periodic. Namely, let us formally calculate the Fourier coef- 


ficients of 


dtr 




V ’ d£ 


U V 1 -v r 

^ J 


2L+ 1 

= _L_ y L, 

2L+] ^ l v’ 9£ f 


/= 1 


dll!, \ -in 1X^3 


(n = —L, ..., L) 


(63) 


Then, it is well-known (see ref. 24 by Kolmogorov and Fomin) that the time-periodic function 


{<%- X , 

n = -L 

minimizes the following functional 




u « 

u v , ^ 


in /t 


2k 


(1 = 0, 1, .... 2L+1) 


(64) 


9u,,l f dv,,^ 




U V’ 9^ 


'v» 
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is a usual Euclidean norm) on the class of periodic functions; i.e., 


r v> 


dv v l 


from formula (64) is 


the best periodic least squares approximation of 


3*0 

u v’-^ 


. Relying on this property, we will further use 


Fourier coefficients (63) as the boundary data that “drive” the ABC’s (which may be referred to as the 
generalized treatment of the boundary data as being time-periodic). Clearly, as we integrate the Navier- 


Stokes equations in time and approach the true oscillatory regime, the “source” function 


u. 


3u v 


and 


its Fourier series 


3v „ ) 


* K 


( :n 3Uv ' 


u v> 3X0 specified on the curve T, which is positioned arbitrarily with respect to 
d ^> ) 


(eq. (64)) also approach each other. 

We now implement the DPM (refs. 9 and 10) to actually calculate the ABC’s. We note that the 
boundary data 

coordinate lines of the grid 9{°. (See formula (53).) Moreover, we do not impose any restrictions on the 
shape of T itself. In our opinion, the DPM (refs. 9 and 10) provides an ideal tool for treating such geo- 
metrically complicated problems. 

Let us introduce the following discrete sets. We consider a six-node two-dimensional stencil 

St m +1/2, j ~ 1 (*m’ + 1 )> ( x m < >j - 1 )»<**+ 1’ V’ ( x m+ 1> Vj + 1 >’ + 1- - 1 > > 

This stencil is actually a projection of the one from figure 6 onto the plane t = const. Obviously, the dis- 
cretization (52) was obtained using St m + 1/2 j. Then, we define 

*t in = 5W° n D- <M = *1= KJ St n 


(x m 


U 

1 / 2 * 


+ 1 /2, j 


°^in ~ kJ St m+\/2,j 

(*m+W2 M,* 


y = *L r\ *Lu 


Clearly, the set of grid nodes y is located near the artificial boundary T. We will call this set the grid 
boundary (by analogy). The sets 9rf in , 9^ in> and y are shown in figure 7. 

Further, we will need to interpolate grid functions from to the points v t <z . Let us select all 
those nodes k c 5\£° that should be taken into account once constructing local interpolation formulas of 
sufficiently high (e.g., second) order. All the nodes k are located not far from Tj. Without loss of 
generality, we always may assume that k c= ry\ We denote the operation of local interpolation from the 
Cartesian grid Vj by R 

Let us also introduce the set of collocation points ocT and the space of eight-component vector 

functions 3 defined on the set a. The elements of will be used as unknowns for the bound- 
ary operator equation, which will replace the exterior linear difference problem. Henceforth, we will 

du n 


treat w£ as vectors containing the values of u n f 0”, p n , and p n and the values of the derivatives 
3v n db n dd n 

-T- 7 T , -4? r , and at the points a; here, C is the (outward) normal to T. Note that the functions w£ 

dC, dg 


a?’ 


are 


the discrete approximations of 


u 


a u ? 

v, "ac” 


n ~ 1 from section 2.1 . 


j 


Generally, the sizes h x and h y of the grid and the size h a of the one-dimensional collocation grid 
on the curve T are not independent; they should be correlated to each other in a certain way. This 
requirement is a consequence of convergence conditions for the DPM algorithm. (See ref. 9.) The 
theoretical questions concerning the correlation between the sizes of grids and a are delineated in 
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Figure 7. Grid sets. Grid 9 {^ — continuous thin vertical lines; 3W° — continuous thin horizontal lines & dashed thin vertical 
lines; continuous boundary f — thick dark line; — gray boxes; M— gray circles; 5\Gn — white boxes; white circles; 

Y = Tin 9 »(_ in — big white circles & boxes. 

reference 9. As concerns practical applications, the final choice of grids is always done taking into 
account some previous computational results. In particular, it seems useful to conduct the computations 
(see refs. 1 1 and 12) for the set of collocation points, which is more concentrated at the outflow part of 
the external boundary in the wake region and uniformly spaced at the inflow part of the external bound- 
ary. Moreover, sometimes the relation |o| - |y| 1/2 appears to be proper. At any rate, for each specific 
class of problems (determined both by the geometry of computational domain and by the parameters of 
fluid at infinity) one always can make an appropriate choice of the grids 5\£° and a relying on general 
theory (ref. 9) and on the numerical experience. 

Let us now specify some wj e W£ and implement the following procedure. First, we smoothly 
interpolate w£ along T (i.e., along the smooth components of O and get the function Rw£ . Here, R is 
an interpolation operator. Then, we drop the normals from the points y to Y and find the values of Rw£ 
at the foots of these normals. Since w£ (and consequently Rw£ ) contains the values of both 0 n ,p n , 

and p" and their normal derivatives and since the distance between any node y and the curve T is small 
(of order h), we may approximately find u n , v n , p'\ and p" at the nodes y using the two first terms of 
the Taylor expansion. We will designate the entire operation of continuation of the boundary data from 
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a to y as 7t wg = u y * Note that the preceding algorithm of continuation generally applies to the 
smooth parts of r (where the normal exists). In practice, however, the curve T is usually not smooth (see 
fig. 1), and it is impossible to construct an appropriate normal when the node y is located in some neigh- 
borhood of the breaking point of the curve. The construction of the operator tl ^ in this case is based on 
the existence of two linearly independent directions along the curve, which enables us to obtain the 
desirable continuation anyway. 

Now, using the calculated continuation of the boundary data, Uy = 7ty a wg , we construct the fol- 
lowing grid function: 


U ^° 


m,; 


u y| if yJ e y 

' l m, j J 

0 if ( x m , yj) e 9i°\y 


(65) 


which is defined already on the entire grid (See formula (53).) Then, we substitute the function u'^-o 
from formula (65) into the left-hand side of system (52). Generally speaking, u^o does not satisfy equa- 
tions (52). Therefore, we generate some nonzero right-hand side, which we designate Lgu^o. Here, Lg 
is the linear operator defined by the left-hand side of system (52). This operator takes the functions 
defined on the grid Of 0 (eq. (53)) as an input and generates the functions defined on the grid Orf 0 
(eq. (54)) as a result. Finally, we truncate the function Lgu^o to the set 0d in , which yields 


f n 

1 n 




m + 1 / 2 , 


Lgu^o 


m + 1 / 2 , j 


0 


if (*m + l/ 2 *J'y) e 
if ( x m + l/2<yj)* M in 


( 66 ) 


We will use the function f^o = fg, + 1 / 2 , j from formula (66) as the right-hand side for the difference AP\ 
by definition, suppf£ + 1 / 2 , j c D in . Once we solve the difference AP with the right-hand side f^o (eq. 
(66)), we get the function Ggf^o. Here, Gg is the Green (i.e., inverse) operator of the difference AP. 
The function Ggf^o is defined on the grid Of 0 . As it is considered only on the sub-grid Of c of 0 , it is 
called the difference potential with the density u” , PJj^Uy = Ggf^oj . (See ref. 9.) Clearly, the differ- 
ence potential satisfies equations (52) since f^o = 0 on 0d\ moreover, it satisfies the boundary condi- 
tions of the difference AP. The difference potential PJ^Uy is a discrete realization of the generalized 

potential mentioned in the introduction. Later, we will find an approximate (i.e., difference) solution to 
the problem (eqs. (8)-(10)) in the form of a difference potential and then use this solution to construct 
the ABC’s, i.e., to obtain the missing relations between the unknowns at vcT and at Vj <z Tj . There- 
fore, we will need to know the difference potential only on the two subsets of Of located closely to T 
and Tj on y c Of and on k c Of, respectively. 

Indeed, once we calculate the difference potential on y, we can then construct the operator P" as 
the trace of the potential, P^iiy *=^P^Uy| . This operator will be the key element of the boundary 

operator equation of the DPM. Actually, P” is a difference boundary projection operator (ref. 9), which 
substitutes Pp (see section 2.1, eq. (1 1)) in practical computations. 
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Once we calculate the difference potential on k, we can interpolate it and find the trace of the solu- 
tion to the linear problem on Vj, = P”^Uy = R v ^P^,Uy . Thereby, we obtain the desirable rela- 
tions between the unknowns at T and F j through the solution of the linearized exterior problem. Let us 
now recall that we have replaced the original infinite-in-space problem by the periodic problem formu- 
lated on the strip {-°° < x < « } x j~ < y < ^ j , claiming that we will need to know the solution only on 

some neighborhood of D in , and therefore the convergence on a fixed interval, -y<y<y, is sufficient 
for our purposes (see subsection 2.3). Since we do not need to calculate the difference potential any- 
where except at y and k, the previous statement is now justified. 

We now formulate the main result of the DPM theory. (See ref. 9.) Consider the entire space of grid 
functions Uy defined on y. Those and only those elements of this space which satisfy the equation 


P"u" = u" 


(67) 


can be complemented to 5\£so that the complement solves system (52) (with boundary conditions (55) 
and (59)). As previously mentioned, the operator P” is a projection; it is a discrete analogue of the Cal- 


deron boundary pseudodifferential operators. (See ref. 17.) For any Uy, the result P"iiy 
P” y Uy) can be explicitly calculated in accordance with the aforementioned 


lit 


’ 


.Lgu^o: 


=> f»0 => G"f: 


n r n 


(TiM 0 


^ r y u Y> r v,Y u r • 


(as well as 
procedure, 


In section 2.1, we have declared our goal as to characterize those and only those 


vp, 


3vf 


that 


would solve equations (8) with boundary conditions (9) on D ex and coincide with the trace of the solu- 
tion on T. Instead, we have provided an analogous classification (see eq. (67)) for the discrete rather 
than for the continuous formulation of the problem. Therefore, we have equivalently replaced linear 
system (52) on along with the boundary conditions (55) and (59), by the boundary equation with 
projection (eq. (67)). Consequently, we can now specify the proper boundary data (see equality (10)) for 


the discrete counterpart of the problem (eqs. (8)— (10)). Namely, let us take any 


du" 


(provided 


from inside D in ) and interpolate it along T to the set of collocation points o, < = R av 

Then, continue w£ using the operator and finally apply P" (which requires solving the AP). In 
accordance with the previously formulated main result, the grid function 


dUy 

,v> aT 


v y = p ;v R ° v 


‘~n 

“v, 


(68) 


admits the complement to 9i that solves the problem (eqs. (52), (55), and (59)). 

Thus, we have completed the first stage of constructing the ABC’s (section 2.1) and now proceed to 
the second one. Instead of the problem (eqs. (8)— (10)), we will consider its discrete counterpart: to solve 
system (52) on 9i with external boundai 7 conditions (55) and (59), and with boundary condition 
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(69) 


v 


n 

y 


at y, Vy in equality (69) comes from formula (68). The solvability of the problem (eqs. (52), (55), (59), 
and (69)) is guaranteed by the special type of boundary data provided in formula (68). 

To actually find the solution to the problem (eqs. (52), (55), (59), and (69)), we calculate the differ- 
ence potential P^. v" with the density v y from formula (68). Then, we interpolate the potential from 
to Vj, which yields 


_ d p n D 

Uv l <*V 


Uv, 


K 


dg jn 


35 J 


(70) 


Equality (70) provides the missing relations between the unknowns at v and at Vj in the Fourier space; 
these relations are based on the solution to the linearized exterior problem. Therefore, equality (70) is 
almost the desirable ABC, and the only remaining step is the inverse Fourier transform in time. Before 
implementing the inverse transform, let us note that the entire algorithm becomes most convenient from 

a practical standpoint if we calculate the matrix representation of the operator T" from formula (70). To 

do that, we choose some basis in W£ , e.g., the simplest one, composed of the vectors like (0, . . . , 0, 1 , 0, 
...,0), and implement the entire proceeding procedure. More precisely, we calculate 

Uv, = R v a wS , for each basis vector w£ . In so doing, we obtain the matrix of 

Ry P ^yy P y ^ yo ( eac h column will be the response to a specific basic function w£ ) and then, multiply- 
ing the above matrix from the right by the interpolation matrix R^, we finally obtain the matrix repre- 


sentation of T" . (Note that we do not start from basis functions 


Uv 


aid"' 

' as 


since the number of nodes c 


is usually much less than the number of nodes v.) In fact, it is possible to show (see ref. 9) that for any 
w a , P^P^Jt yo wS = PlyyKyowS and therefore we need to calculate the matrix Ry^P^^g • Clearly, 


the computation of each column of the matrix R v ^P^n yo requires one solution of the difference AP, 

which, in turn, involves the direct (eq. (56b)) and inverse (eq. (61)) Fourier transforms and the solution 
of the problem (eqs. (57) and (59)) for each wavenumber k,k = -J, Either Fourier transform will 

require only 0(M J) rather than 0(M J 2 ) operations. (For definitions of M and J, see formula (46).) 
Indeed, the support of the right-hand side f^o is actually concentrated near V since u^o differs from 
zero only on y and the operator Lg is local. Therefore, while calculating direct Fourier transform (56b) 

for each m, m = 0, 1 M - 1 , only a few values C + 1 / 2 , j differ from zero, and consequently the total 

cost of this computation is 0(M ■ J) operations. Analogously, while calculating the inverse Fourier 
transform (61) for each m, m = 0, 1, ..., M, we need to know uJJ I> j only for a few selected values of j 
since all other (x m , yj) do not belong to K. Therefore, the total cost of this computation is 0(M J) 
operations as well. Finally, the solution of the problem (eqs. (57) and (59)) for each wavenumber k , 
k — —J, ...,J costs 0(M) operations. (See ref. 24.) Adding all these quantities, we obtain a total of 

0(M J) operations for the computation of each column of the matrix R v ^P^,7t ya . We see that 

although the entire algorithm requires repeated solution of the difference AP, the solution may be 
obtained by means of an efficient procedure, which should make the total expense for calculating the 
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ABC’s quite acceptable. Note that in our previous work (see refs. 8, 1 1, and 12) we have used a differ- 
ent version of the algorithm. We expect that the total cost per one Fourier mode in time will decrease for 
the current version because of using the thin-layer rather than the full Navier-Stokes equations. (Indeed, 

the matrices A£ and B£ in system (57) are 4 x 4 and the matrices for the full Navier-Stokes equations 
are 8 x 8.) (See ref. 8.) 


Recall, our final goal is to express the values of physical variables at Vj, u y s ( u v , v , p , p ), 
( 3u„ > ] 

• Choosing the same discretization in time as in the formula (eq. (62)) and imple- 


in terms of 


u v’ ^ 


menting inverse Fourier transform (64), we obtain from formula (70), 


n - L 


j/i 


n = -L 


du^ 

Kf 


in lx 


2k 


(/ = 0, 1, ..., 2L + 1) 


(71) 


Then, substituting expression (63) into formula (71) and changing the summation order, we get, 


n = L . 2jt 2Z.+ 1 / -n s 

= y f n e T_L_ y L ^ 

v i ^ 2L + 1 ^ v> ac 

n ~ -L s = 1 V 
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L 
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Finally, designating 


jl, S — 


, n ~ L . N 2k \ 

1 ^ 


2 L + 


7 I 

n - -L 


(1 = 0, 1, ..., 2L + 1) 


we obtain 


2L+i ( du^ 


<- X T '' 


5 = 1 


II •> 

V v a? , 


(/ = 0, 1, ..., 2L + 1) 


(72) 


Equality (72), which is a specification of equality (11), provides the missing boundary relations between 
the values of the unknowns at T r and at rf (in the discrete formulation). Therefore, equality (72) pro- 
vides the ABC’s we were aiming to obtain. Additionally, we note that the ABC’s (eq. (72)) can be sim- 
plified for the case of integrating the Navier-Stokes equations step-by-step in time inside D in . In doing 
so, we only need to know u v on the upper time level, i.e., for t=T, which corresponds to / = 2L + 1 or 
to / = 0 because of the periodicity. Substituting / = 0 into equality (72), we obtain 


2L + ] 


I t = T 


- X T °' 


i = 1 


3u* 


(73) 


Equality (73) is a desirable global ABC for implementation together with the step-by-step integration 
procedure in time. Indeed, formula (73) expresses the values of u, v, p, and p (perturbations) at the 
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outermost coordinate row v | on the upper time level t = T as a function of the prescribed data 


du v 


through the time-periodic solution of the linearized thin-layer equations with the free-stream boundary 
condition at infinity. We note that the matrices of operators T 0, 5 are calculated explicitly and therefore, 
the practical implementation of the ABC’s (eq. (73)) is reduced to a few matrix-vector multiplications. 
We also note that this practical implementation may preliminarily require some interpolation in time at 
the nodes v, which may be represented in a matrix form as well. If we use an explicit scheme for inte- 
grating the Navier-Stokes equations inside D in , then we directly implement formula (73) at each time 
step for determining the missing values of the unknowns at the outermost coordinate row of the grid on 
the upper time level. If the scheme inside D in is implicit, then we include the relations (eq. (73)) into the 


system of equations we solve on the upper time level treating all T 0, 


U“ 


ac 


for s < 2L + 1 as forcing 


terms. 


4. Concluding Remarks 

We have constructed the DPM-based nonlocal artificial boundary conditions for computation of 
oscillating external flows, specifically, compressible viscous fluid flows past finite bodies. To develop 
the ABC’s, we used linearization of the governing equations against the constant free-stream back- 
ground in the far field. To justify the constructions of difference potentials used for computation of the 
ABC’s, we provided some results on solvability of the linearized thin-layer equations. The nonlocal 
nature of the proposed ABC’s arises from their closeness to the exact boundary conditions. In spite of 
this nonlocal nature, our ABC’s apply to artificial boundaries of irregular shape with equal ease, which 
is very important for applications. We expect that these boundary conditions may become an effective 
numerical tool in practical computing. The numerical results on the implementation of the described 
ABC’s will be presented in future work. 

Note that we describe the algorithm for calculating the ABC’s only for a particular class of methods 
used for integrating the Navier-Stokes equations inside D in , namely, for such methods that the knowl- 
edge of missing relations between only two external coordinate rows of the grid (v and V]) is sufficient 
for closing the discrete system inside the computational domain. Obviously, once the method used 
inside D in is of higher (than the second) order, the consideration of only two curves, T and T ] , might be 
insufficient. However, we always can assume that the “linear region” D ex contains more than one curve, 
e.g., T} and T 2 instead of only T^ and can treat this case in the same way as described in this paper. 
Moreover, one can use higher order schemes for solving the linearized exterior problem as well. Such 
modifications may extend the possible range of applications for the described technique by including, 
for example, some computational problems of aeroacoustics. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
April 11, 1996 
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